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INTRODUCTION 


In  recent  years,  studies  in  soil  dynamics  have  dealt  extensively 
with  wave  propagation  in  soils.  Both  laboratory  and  field  wave  veloci¬ 
ties  have  been  measured  in  an  effort  to  determine  the  various 
parameters  that  affect  wave  propagation.  Wave  propagation  is  an  ex¬ 
cellent  technique  for  applying  extremely  small  strains  to  soils  and  to 
accurately  measure  the  stress-strain  response  of  the  soils  due  to  the 
small  strains.  Such  small  strains  are  nearly  elastic.  Consequently, 
studies  of  wave  propagation  in  soil  'dynamics  have  contributed  greatly  to 
the  understanding  of  elastic  soil  behavior. 

However,  due  to  the  use  of  limiting  equilibrium  analyses  in  soil 
statics,  studies  in  soil  statics  have  tended  to  concentrate  on  the 
strength  of  soils  and  on  large  strain  behavior.  Large  strains  involve 
slippage,  rearrangement,  and  crushing  of  particles.  These  are  sources 
of  plastic  strain.  The  studies  of  plastic  strains  in  soil  statics  have 
led  to  the  use  of  work  hardening  plasticity  theories  in  soil  mechanics, 
the  development  of  critical  state  soil  mechanics,  and  to  the  development 
of  stress-dilatancy  theory  (Drucker,  et  al.  1957;  Rowe,  1962;  and 
Schofield  and  Wroth,  1968). 

Hardin  (1978)  presented  a  state-of-the-art  paper  on  the  stress- 
strain  behavior  of  soils  in  which  information  on  the  elastic  behavior  of 
soils  from  soil  dynamics  is  synthesized  with  information  on  the 
plastic  behavior  of  soils  from  soil  statics.  Synthesizing  these  two 
types  of  soil  behavior  provided  a  framework  for  the  development  of 
comprehensive  elasto-plastic  stress-strain  relationships  for  soils.  The 
stress-strain  relations  are  formulated  in  terms  of  effective  stress  and 
particulate  mechanics  is  used  as  a  guide  for  the  development  of 
macroscopic  constitutive  equations.  In  particular,  the  Hertz  and 
Mindlin  (1949)  theories  have  to  do  with  the  elastic  behavior  of  soils 
and  the  stress-dilatancy  theory  is  related  to  the  plastic  behavior  of 
soils  (Rowe,  1971). 

Many  constitutive  equations  currently  being  proposed  for  soils 
arise  from  studies  in  continuum  mechanics.  They  sure  often  csurefully 
constructed  from  the  continuum  mechanics  point  of  view,  but  show  little 
regard  for  the  constraints  imposed  by  particulate  mechanics,  i.e.  the 
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kinematics  of  particle  movement,  particle  crushing,  and  the  bonding  at 
particle  contacts.  As  a  result,  the  physical  meaning  of  the  coeffi¬ 
cients  in  the  constitutive  equations  is  often  obscure,  and  a  given  set 
of  values  for  the  coefficients  is  valid  for  a  very  restricted  range  of 
loading. 

To  provide  physically  based  elasto-plastic  soil  constitutive 
equations,  the  U.  S.  Air  Force  Office  of  Scientific  Research  has 
sponsored  the  research  summarized  in  this  report  through  its  Project  No. 
AP0SR-84-0195.  The  overall  project  objective  is  the  development  and 
verification  of  three-dimensional  elasto-plastic  constitutive  equations 
for  soils.  Essential  features  of  soil  behavior  that  result  from  the 
soil  skeleton  being  particulate  are  included  and  the  developed  soil 
model  incorporates  the  oehavior  characteristics  indicated  by  particulate 
mechanics  theories.  The  construction  of  the  elastic  component  of  the 
constitutive  equations  is  influenced  by  the  Hertz  and  Mindlin  theories. 
The  plastic  component  is  formulated  using  the  stress-dilatancy  theory 
and  includes  the  dependency  of  particulate  materials  on  the  effective 
stress  increment  direction  as  well  as  state  of  effective  stress.  This 
is  accounted  for  by  defining  two  plasticity  models  for  two  respective 
classes  of  stress  increment  directions  with  different  plastic  potential 
and  hardening  functions  for  each  class. 

Because  the  developed  constitutive  equations  have  been  constructed 
by  synthesizing  small  strain  data  from  soil  dynamics  with  data  for  large 
strains  from  soil  statics  (Hardin,  1978),  a  given  set  of  coefficients  is 
valid  for  a  wide  range  of  loading  conditions,  from  wave  propagation 
strain  levels  to  loadings  where  the  soil  fails.  The  constitutive  equa¬ 
tions  include  such  state  and  stress  history  parameters  as  void  ratio  and 
overconsolidation  ratio.  Consequently,  some  of  the  effects  of  changes 
in  these  parameters  are  accounted  for  directly  and  therefore  do  not 
cause  variations  in  the  coefficients.  Many  of  the  coefficients  are 
nearly  constant  for  a  wide  variety  of  soils  and  loading  conditions,  and 
the  effects  of  soil  disturbance  on  some  of  the  coefficients  is  beginning 
to  be  understood. 

This  report  focuses  on  the  developed  elasto-plastic  constitutive 
equations  and  their  implementation  into  a  three-dimensional  finite 
element  code  EPSAP  (Elasto-Plastic  Soil  Analysis  Program) .  The  deve- 
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loped  constitutive  equations  are  based  on  assuming  the  soil  to  be  a 
continuum.  This  is  necessary  in  order  to  obtain  solutions  to  practical 
geotechnical  engineering  problems.  However,  one  should  always  remember 
that  the  soil  is  composed  of  discontinuous  particles  and  the  properties 
of  the  corunuun  should  reflect  the  essential  features  of  soil  behavior 

A 

that  result  from  the  particulate  nature  of  soils.  Primary  emphasis  of 
this  project  has  been  on  the  understanding  and  realistic  modeling  of  the 
constitutive  behavior  of  soils.  The  finite  element  program  has  been 
developed  for  use  in  the  verification  of  the  constitutive  model;  but  can 
also  be  used  for  the  analysis  of  three-dimensional  soil  systems. 

At  the  beginning  of  the  project  in  1984,  the  principal  investi¬ 
gators  believed  that  the  constitutive  equations  by  Hardin  (1978)  could 
be  perfected  and  implemented  into  a  three-dimensional  finite  element 
analysis  program  in  a  relatively  short  time.  It  would  then  be  possible 
to  solve  boundary  value  problems  corresponding  to  laboratory  tests,  e.g. 
the  vertical  loading  of  model  footings  on  sands.  However,  the  ideas  for 
improvement  of  constitutive  models  that  have  been  discovered  during  the 
project  were  unforseen  and  the  complexity  of  numerical  implementation 
was  underestimated.  Consequently,  as  the  project  concludes,  the  ana¬ 
lysis  program  is  just  now  ready  for  solution  of  problems  involving 
monotonic  loading  of  cohesionless  soils.  The  development  of  consti¬ 
tutive  equations  for  cohesive  soils  and  cyclic  loading  has  been  advanced 
but  requires  considerable  additional  effort.  Thus,  the  verification  of 
constitutive  equations  has  been  by  comparison  to  test  results  for 
specimens  subjected  to  homogenous  states  of  stress  and  strain.  Veri¬ 
fication  by  finite  element  simulation  of  footing  tests,  etc.,  has  not 
been  completed. 

A  complete  set  of  constitutive  equations  and  numerical  procedure 
for  the  monotonic  load  analysis  of  cohesionless  soils  are  presented 
herein  and  will  be  submitted  for  publication  in  the  near  future.  Though 
the  focus  is  on  the  final  stage  of  the  supported  research  effort,  some 
of  the  other  accomplishments  previously  reported  are  summarized.  A 
brief  overview  and  summary  of  the  published  results  for  the  project  are 
presented  in  the  next  section. 
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PROJECT  OVERVIEW  AND  PUBLICATION  SUMMARIES 

Three-dimensional  elasto-plastic  constitutive  equations  for  soils 
are  comprised  of  the  elemental  parts  listed  in  Table  1.  Much  of  the 
framework  in  Table  1  and  preliminary  definition  of  many  of  the  elemental 
parts  was  presented  by  Hardin  (1978).  Research  completed  under  the 
current  contract  has  aimed  primarily  toward  improved  definition  of  the 
plastic  potential  surfaces  and  work  hardening  and  softening  functions 
for  monotonic  loading.  An  important  development  discovered  during  the 
past  year  is  the  implicit  definition  of  plastic  potential  surfaces. 
Principal  plastic  strain  increment  ratio  functions  representing  the 
gradients  of  plastic  potential  functions  are  defined.  This  provides 
flexibility  in  defining  surface  shapes  that  accurately  model  soil 
behavior  and  simplifies  numerical  evaluation  of  the  plastic  strain 
increments . 

Details  of  accomplishments  have  been  reported  to  APOSR  in  project 
reports  and  through  copies  of  papers  published  or  submitted  for  publi¬ 
cation.  The  items  studied  are  listed  in  Table  2  with  references  to 
appropriate  papers  or  reports  containing  resultB.  A  brief  summary  of 
each  paper  submitted  for  publication  is  included  below  in  chronological 
order. 

Soil  Particle  Crushing  (Hardin,  1985a).  -  The  strength  and  stress-strain 
behavior  of  an  element  of  soil  is  affected  greatly  by  the  degree  to 
which  particle  crushing  or  particle  breakage  takes  place  during  loading 
and  deformation.  For  the  type  of  deformation  that  primarily  produces 
volume  change,  such  as  one-dimensional  strain  or  isotropic  compression 
(Class  2  stress  increments),  particle  crushing  adds  to  the  reduction  in 
volume.  For  the  type  of  deformation  where  particles  sue  moving  past  or 
around  one  another,  such  as  triaxial  compression  or  simple  shear  (Class 
1  stress  increments),  crushing  at  sliding  contacts  decreases  the  rate  of 
dilation  corresponding  to  a  given  principal  stress  ratio. 

In  order  to  understand  the  physics  of  the  strength  and  stress- 
strain  behavior  of  soils  and  to  devise  mathematical  models  that  ade¬ 
quately  represent  such  behavior,  it  is  important  to  define  the  degree  to 
which  the  particles  of  an  element  of  soil  are  crushed.  The  amount  of 
particle  crushing  that  occurs  in  an  element  of  soil  under  stress  depends 


4 


EFFECTS  OF  SOIL  DISTURBANCE  AND  AGING 


I 


Table  1.-  ELEMENTAL  PARTS  OF 
3D  CONSTITUTIVE  EQUATIONS 


*  ELASTIC  COMPONENT  (Hertzian) 


*  PLASTIC  COMPONENT 

:  DEFINITION  OF  LOADING 

;  For  Monoton 

Class  1 

PLASTiC  POTENTIAL 
SURFACE 

ic  Loading 

Class  2 

PLASTIC  POTENTIAL 
SURFACE 

:  WORK  HARDENING  AND 
:  SOFTENING  FUNCTION 

WORK  HARDENING 
FUNCTION 

For  Cyclic 

Class  1 

YIELD  RATE 
:  SURFACE 

Loading 

Class  2 

PRECONSOLIDATION 

SURFACE 

V 

MASING 

V 

V 

UNLOAD-RELOAD 

SHAPE 

;  KINEMATIC  HARDENING 

KINEMATIC  HARDENING 

:  ANISOTROPY  (PACKING  AND  STRESS) 

:  ROTATION  OF  STRESS  INCREMENT 

> 

*  VISCOUS  COMPONENT  (strain  Rate  Effects) 
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Table  2  -  Research  Accomplishments 


Constitutive  Behavior 


Task 

Soil  Elasticity 


Particle  Crushing 


Papers.  Reports  and  Comments 

Hardin,  B.  0.  and  Blandford,  G.  E.  (1988a), 
"Elasticity  of  Particulate  Materials," 
submitted  to  Journal  of  Geotechnical 
Engineering:  draft  submitted  to  APOSR. 

Hardin,  B.  0.  (1985),  "Crushing  of  Soil 
Particles,"  Journal  of  Geotechnical 
Engineering,  ASCE,  Vol.  Ill,  No.  10,  pp. 
1177-1192. 


Strength  of  Soils  Hardin,  B.  0.  (1985),  "Strength  of  Soils  in 

Terms  of  Effective  Stress,"  Proceedings  of 
the  Richart  Commemorative  Lectures,  ASCE, 
Detroit,  October,  pp.  1-78. 

Hardin,  B.  0.  (1988b),  "Tbe  Low-Stress 
Dilation  Test,"  accepted  for  publication  in 
the  Journal  of  Geotechnical  Engineering, 
ASCE:  draft  submitted  to  APOSR. 

Class  1  Plastic  Potential  "Class  1  Plastic  Potential  Function  Formu- 
Surface  lation  Appendix  III,  first  annual  report 

to  AFOSR,  Sep.  26,  1985. 


Class  2  Plastic  Potential:  Modification  of 
the  Class  1  plastic  Surface  potential 
function  also  applies  to  Class  2, because 
the  Class  1  function  is  a  component  of  the 
Class  2  function. 


Class  1  Work  Hardening  "Modeling  Work  Softening  Behavior  for  Class 

and  Softening  Function  1  Plastic  Hardening,"  Appendix  II,  second 

annual  report  to  APOSR,  Sep.  15,  1986. 

Hardin  (1988d),  "Triaxial  Compression  for 
Cohesionless  Soils,"  to  be  submitted  to 
Journal  of  Geotechnical  Engineering,  ASCE. 
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Table  2  -  Research  Accomplishments  (continued) 


Task 

Class  1  Work  Hardening 
and  Softening  Function 
(continued) 

Class  2  Work  Hardening 
Function 


Application 


Papers,  Reports  and  Comments 

A  database  which  contains  several  hundred 
stress-strain  curves  from  drained  triaxial 
compression  testa  has  been  completed. 

"Modeling  Class  2  Plastic  Hardening," 
Appendix  III,  second  annual  report  to  APOSR, 
Sep.  15,  1986. 

Hardin  (1987),  "ID  Strain  in  Normally 
Consolidated  Cohesionless  Soils,"  Journal 
of  Geotechnical  Engineering,  ASCE,  Vol.  113, 
No.  12,  pp.  1449  -  1467. 

Hardin  (1988a),  "ID  Strain  in  Normally 
Consolidated  Cohesive  Soils,"  Journal  of 
Geotechnical  Engineering,  ASCE,  39  pp. ,  to 
be  published. 

Hardin,  et  al.  (1988),  "Triaxial  Compres¬ 
sion,  Simple  Shear,  and  Strength  of  Wheat, 
International  Sumner  Meeting  of  the  ASAE, 
Rapid  City,  SD,  June  26-28,  1988,  Paper 
#88-4022;  in  review  for  possible  publica¬ 
tion  in  the  Transactions  of  the  ASAE. 


Finite  Element  Analysis 


Task 


Papers,  Reports  and  Comments 


Matrix  Formulation  of 
Constitutive  Equations 


"Matrix  Formulation  of  Constitutive 
Equations,"  section  of  first  annual  report 
to  APOSR,  Sep.  26,  1985,  pp.  8-13. 


3D  Constitutive  Equations  Hardin  and  Blandford  (1988b),  "3D  Consti¬ 
tutive  Equations  for  Cohesionless  Soils," 
to  be  submitted  to  Journal  of  Geotechnical 
Engineering,  ASCE. 


Implementation  of  Because  the  constitutive  equations  are 

Constitutive  Equations  being  developed  and  revised,  the  finite 

into  Finite  Element  Code  element  code  must  be  continuously  updated 

to  include  newly  developed  models. 
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Table  2  -  Research  Accomplishments  (continued) 

Task 

Pacers,  Reports  and  Comments 

Development  of  Three- 
Dimensional  Element 

Library 

"Finite  Element  Formulation  Details," 
Appendix  V,  second  annual  report  to  AFOSR, 
Sep.  15,  1986.  (includes  representation 
of  infinite  boundaries);  plus  this  report. 

Investigation  of  Nonlinear 
Solution  Strategies 

"Nonlinear  Solution  Algorithm,"  section  of 
second  annual  report  to  AFOSR,  Sep.  15, 

1986,  pp.  16-22;  and  this  report. 

Finite  Element  Code 
Verification 

A  soil  specimen  subjected  to  homogeneous 
states  of  stress  along  a  specified  stress 
path  is  represented  using  one  modified 
hexahedron  finite  element  for  the  symmetric 
one-eighth  of  the  problem.  The  computed 
stress-strain  relation  is  compared  to  that 
computed  directly  from  the  constitutive 
equations . 

Some  results  of  the  finite  element  part  of 
this  verification  procedure  are  given  in 
"Verification  Problems,"  section  of  second 
annual  report  to  AFOSR,  Sep.  15,  1986, 
pp.  22-27;  and  this  report. 

3D  Elasto-Plastic  Finite 
Element  Analysis 

Blandford  and  Hardin  (1988),  "3D  Elasto- 
Plastic  Finite  Element  Analysis  of 
Cohesionless  Soils,"  to  be  submitted  to 
International  Journal  for  Numerical  and 
Analytical  Methods  in  Geomechanics. 
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on  particle  size  distribution,  particle  shape,  effective  stress  state, 
effective  stress  path,  void  ratio,  particle  hardness  (i.e.  hardness  of 
cementing  material  or  weakest  constituent  of  a  particle  and  weakest 
particles  of  an  element),  and  the  presence  or  absence  of  water. 

Data  for  a  large  nimber  of  single  mineral  soils  and  rockfill-like 
materials  have  been  analyzed  and  equations  developed  that  can  be  used  to 
estimate  the  total  breakage  expected  for  a  given  soil  subjected  to  a 
given  loading.  In  doing  this,  new  measures  of  particle  breakage  have 
been  defined,  called  breakage  potential,  total  breakage,  and  relative 
breakage,  that  integrate  the  changes  in  the  particle  size  distribution 
curve  for  all  sizes  greater  than  0.074  mm.  A  methodology  based  on  pro¬ 
bability  analysis  for  studying  the  breakage  of  particles  of  each  size 
within  a  distribution  of  sizes  has  been  outlined.  Continued  development 
of  the  probability  model  for  particle  crushing  is  currently  supported  by 
the  National  Science  Foundation. 

Soil  Strength  in  Terms  of  Effective  Stress  (Hardin,  1985b).  -  An  under¬ 
standing  of  the  strength  of  soils  is  essential  to  the  formulation  of 
general-purpose  constitutive  models  for  soils.  Failure  is  a  limiting 
case  and  should  be  defined  first,  so  that  the  form  of  equations  adopted 
for  various  elements  of  the  constitutive  model  (the  Class  1  hardening 
function  in  particular)  will  behave  properly  as  limiting  conditions  are 
approached.  The  Mohr-Coulomb  failure  theory  is  used  by  most  geotechical 
engineers  to  model  the  strength  of  soils.  However,  a  more  fundamental 
approach  that  considers  the  particulate  nature  of  soils  has  been  pur¬ 
sued. 

Factors  that  determine  soil  strength  (e.g.  confinement,  density, 
particle  size  distribution,  mineral  friction,  aging,  particle  crushing, 
etc.)  produce  their  effect  through  their  influence  on  two  fundamental 
mechanisms  of  resistance  to  deformation  in  particulate  materials.  These 
two  fundamental  mechanisms  are 

*  bonding  at  particle  contacts,  and 

*  the  kinematics  of  particle  movement  within  an 
element  of  deforming  soil. 
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The  first  of  these  fundamental  mechanisms  is  physio-chemical  in  nature 

and  the  second  is  geometric.  Inclusion  of  these  two  basic  mechanisms 

must  be  considered  in  determining  the  appropriate  soil  strength  model. 

Figure  1  illustrates  the  developed  strength  model.  The  two  parameters 

that  represent  the  two  basic  mechanisms  of  soil  strength  are  the  bonding 

obliquity  angle  <p  and  the  maximun  rate  of  dilation  d  .  Variation 
o  max 

of  these  parameters  with  minor  principal  stress  is  illustrated  in 

Figs.  1(b)  and  1(c)  and  the  resulting  strength  envelope  is  shown  on  the 
Mohr  diagram  in  Fig.  1(a).  The  maximum  rate  of  dilation  curve  is  de¬ 
fined  directly  from  the  results  of  strength  tests  (normally  triaxial 
compression) ;  whereas  the  strength  model  must  be  used  to  determine  the 
bonding  obliquity  angle  from  the  same  tests.  Curvature  of  the  strength 
envelope  in  Fig.  1(a)  results  from  variation  with  stress  of  the  two 

fundamental  parameters  4>  and  d 

^  o  max 

The  bonding  obliquity  angle  is  equal  to  the  mineral  friction  angle 
for  cohesionless  soils,  but  includes  cohesive  and  frictional  bonding  for 
cohesive  materials.  A  method  for  isolating  the  magnitudes  of  cohesive 
and  frictional  bonding  in  cohesive  materials  is  proposed.  Analysis  of 
triaxial  compression  test  results  for  stabilized  sands  shows  that  the 
dimensionless  contact  cohesion  parameter  is  a  measure  of  the  strength  of 
cementation  bonding. 

The  model  is  generalized  with  respect  to  stress  path  through  the 
dilatancy  factor  and  the  interpolation  function  F(b)  (refer  to  the  paper 
for  precise  definition  of  these  terms).  The  dilatancy  factor  ties 
together  the  values  of  maximum  rate  of  dilation  for  different  stress 
paths.  Experiments  indicate  that  the  dilatancy  factor  is  independent  of 
stress  path  for  a  given  material  at  a  given  density.  The  interpolation 
function  accounts  for  the  effect  of  stress  path  on  energy  transmission 
and  dissipation  at  the  peak  principal  stress  ratio.  The  interpolation 
function  F(b)  may  be  chosen  to  give  identical  or  differing  strengths  in 
triaxial  compression  and  extension.  It  can  also  be  chosen  so  that  the 
value  of  b  =  ( a £  -  Og)/(oj  -  o'^)  for  plane  strain  matches  test  results. 

Model  calibration  has  been  performed  using  the  results  of  over  700 
triaxial  compression  tests,  most  on  cohesionless  materials.  The 

analyzed  database  includes  silts,  uniform  and  well  graded  sands,  and 
rockfill-like  materials  with  maximun  particle  sizes  up  to  200  am,  as 
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well  as  normally  and  heavily  overconsolidated  clays  and  stabilized 
sands.  Cohesionless  materials  range  in  hardness  from  chlorite  and 
calcite  to  aluninun  oxide,  with  particle  shapes  from  round  to  angular. 
The  minor  principal  stress  varies  from  0.01  to  677  atmospheres,  with  an 
overall  void  ratio  range  from  0.17  to  1.28. 

An  example  is  presented  in  the  paper  showing  how  the  model  may  be 
used  to  estimate  the  strength  of  materials  that  cannot  be  tested.  Thus, 
the  model  may  prove  useful  to  engineers  involved  in  the  design  of  large 
earth  and  rockfill  structures. 

One-Dimensional  Strain  in  Soils  (Hardin,  1987;  1988a).  -  The  first  paper 
compares  a  stress-strain  model  for  one-dimensional  strain  in  normally 
consolidated  cohesionless  soils  to  test  results  covering  the  stress 
range  0  to  1300  atmospheres.  This  model  is  called  the  "1/e  vs.  o^P 
model"  because  the  implied  relationship  between  void  ratio  e  and 
vertical  effective  stress  o^.  is  linear  in  the  low-stress  range.  The 
relationship  between  1/e  and  is  nonlinear  at  higher  stresses  where 
the  effect  of  particle  crushing  becomes  significant.  Effects  of  initial 
void  ratio,  relative  density,  particle  shape,  particle  size 
distribution,  and  particle  material  on  lD-strain  behavior  have  been 
examined  in  the  framework  of  the  model.  The  model  can  accurately  repre¬ 
sent  lD-strain  in  rockf ill-like  materials  as  well  as  sands  and  gravels. 
The  data  and  equations  presented  provide  a  means  of  estimating  shrinkage 
during  placement  of  rockf ills  and  estimating  settlement  of  structures  on 
rockfills  where  it  is  not  practical  to  conduct  lD-strain  tests  or  pene¬ 
tration  tests. 

The  second  paper  compares  a  stress-strain  model  for  one-dimensional 

strain  in  normally  consolidated  cohesive  soils  to  test  results  covering 

the  stress  range  0  to  950  atmospheres.  For  cohesive  soils  the  1/e  vs. 

o;P  relationship  is  linear  in  the  normally  consolidated  range.  The 

slope  of  the  linear  relationship  is  nearly  independent  of  soil 

disturbance  and  its  variation  is  relatively  small  for  a  wide  variety  of 

soils  with  different  compressibilities.  Because  of  this  the  effects  of 

aging  and  soil  disturbance  are  isolated  to  a  single  model  parameter, 

providing  a  quantative  measure  of  these  effects.  The  liquid  limit  is 

easily  related  to  the  model  because  the  a’  =  0  state  is  included.  These 

v 
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features  provide  for  better  interpretation  of  test  results.  An 
undisturbed  sample  can  be  tested  as  usual  and  then  remolded  to  a  water 
content  slightly  above  the  liquid  limit.  The  liquid  limit  should  be 
measured  as  water  is  added  and  the  soil  remolded.  The  final  remolded 
sample  should  then  be  subjected  to  a  ID  strain  test.  Results  for  both 
undisturbed  and  remolded  samples  define  the  model  slope  and  the  liquid 
limit  and  remolded  test  results  are  beneficial  in  choosing  the  best 
value  of  p.  By  following  this  procedure  all  three  tests  are  conducted 
on  identical  material. 

The  one-dimensional  strain  models  for  normally  consolidated  soils 
are  used  to  define  Class  2  plastic  hardening  for  monotonic  loading  (see 
Table  1).  In  order  to  define  Class  2  plastic  hardening  for  cyclic 
loading  the  unload-reload  shape  for  lD-strain  must  be  modeled.  This 
introduces  additional  complexity  because  the  1/e  vs.  o'*5  relationship 
for  unload- reload  is  nonlinear,  the  shape  of  the  unload-reload  curve 
depends  on  stress  history  as  represented  by  the  preconsolidation  stress, 
and  the  preconsolidation  stress  is  just  one  point  on  a  preconsolidation 
surface  when  we  generalize  to  three-dimensions. 

Application  of  Constitutive  Models  to  Wheat  (Hardin,  et  al.,  1988). 
Some  aspects  of  constitutive  behavior  are  common  to  all  particulate 
materials,  e.g.  the  importance  of  bonding  at  particle  contacts  and  the 
dilatancy  of  dense  particulate  materials.  Thus ,  constitutive  equations 
developed  for  cohesionless  soils  can  often  be  applied  to  wheat  with 
relatively  little  modifications  provided  the  material  constants  are 
measured  for  wheat . 

Hardin,  et  al.  (1988)  utilize  models  developed  for  soils  to  define 
the  constitutive  behavior  of  bulk  wheat .  Two  models  define  three- 
dimensional  behavior:  a  model  for  wheat  elasticity  (small-strain 

behavior)  and  a  model  for  the  strength  of  wheat.  An  elasto-plastic 
model  for  the  behavior  of  wheat  in  triaxial  compression  is  presented. 
The  model  includes  the  volume  change  behavior  during  triaxial 
compression  and  is  in  that  sense  three-dimensional.  Simple  shear  and 
triaxial  compression  behaviors  are  compared  by  normalizing  the  two 
moduli  versus  Btrain  relationships. 
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Low  Stress  Dilation  Testing  (Hardin,  1988b;  to  be  published) .  -  The  low- 
stress  dilation  test  is  a  new  test;  used  in  conjunction  with  other  tests 
to  define  the  peak  strength  of  a  particulate  material.  It  provides  a 
measure  of  kinematics  of  particle  movement  within  a  deforming  element  of 
particles  making  its  results  fundamental  to  the  definition  of  soil 
strength  and  to  the  modeling  of  triaxial  compression  in  soil.  The  low- 
stress  dilation  test  relates  to  the  geometric  component  of  soil 
strength . 

The  objective  of  the  test  is  to  measure  the  rate  of  dilation  in  a 

particulate  specimen  deforming  with  minor  principal  stress  =  0. 

Initially,  the  test  was  referred  to  as  an  "unconfined"  dilation  test. 

However,  in  order  to  maintain  necessary  contact  between  specimen  and 

apparatus  the  value  of  <jg  must  be  somewhat  greater  than  zero;  leading  to 

the  name  "low-stress"  dilation  test. 

Dilation  refers  to  the  volume  changes  that  occur  during  deformation 

of  a  particulate  material.  The  relationship  between  volumetric  and 

major  principal  strains  during  triaxial  compression  or  plane  strain 

loading  is  illustrated  schematically  in  Fig.  2.  Rate  of  dilation  d,  as 

measured  by  the  test,  is  the  rate  at  which  volumetric  strain  v  changes 

with  major  principal  strain  (d  =  -dv/d£j).  The  rate  of  dilation  in 

Fig.  2  changes  continuously  with  increasing  €^,  first  increasing  to  the 

maximum  value  d  ,  then  decreasing  to  zero  as  the  strain  becomes  very 
max 

large . 

Maximum  rate  of  dilation  increases  with  decreasing  confinement. 

The  maximum  value  of  d  at  a l  -  0  is  denoted  d  (Fig.  2).  The  low- 

max  3  ocr 

stress  dilation  test  has  been  conceived  as  a  simple  means  of  measuring 
d  . 

oo 

Parameter  d  is  one  of  the  coefficients  in  a  strength  model  pro¬ 
posed  by  Hardin  (1985b),  The  nature  of  the  model  is  illustrated  in  Fig. 
3  where  it  is  applied  to  peak  strength  of  a  cohesionless  rockfill-like 
material  measured  by  four  drained  triaxial  compression  tests.  Two  model 
parameters  represent  the  two  basic  mechanisms  of  soil  strength.  They 
are  the  bonding  obliquity  angle  ♦  ,  which  for  cohesionless  soils  is 
equal  to  the  mineral  friction  angle  ♦  (Fig.  3b),  and  the  maximum  rate 

of  dilation  d  (Fig.  3c).  Both  ♦  and  d  vary  with  oL  as  shown  in 
max  p  max  3 

Fig.  3,  with  maximum  values  occur ing  at  oL  -  0  (♦  =39.8*  and  d  = 

3  pOo  o o 
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Fig.  3.  Strength  Model  and  Test  Results  for  Nap a  Basalt 
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0.71,  respectively).  Hie  curved  Mohr  envelope  defined  by  the  ♦  and 

d  variations  is  compared  to  the  measured  Mohr  circles  in  Fig.  3(a). 
max 

Derivation  of  the  strength  model  and  discussion  of  its  application  to 
cohesive  soils  can  be  found  in  Hardin  (1985b). 

The  test  apparatus  described  in  this  paper  is  relatively  small  in 
size,  designed  to  test  sands.  However,  an  important  feature  of  the  low- 
stress  dilation  test  is  its  potential  for  testing  rockfill-like  materi¬ 
als.  It  should  be  relatively  easy  (compared  to  large-scale  triaxial 

compression  testing)  to  construct  a  low-stress  dilation  apparatus  that 

3 

could  test  1  to  2  id  specimens. 

This  paper  includes  discussion  of  the  strength  model,  description 
of  the  test  apparatus  and  procedure,  and  a  theoretical  analysis  and  test 
results  for  a  packing  of  brass  rods.  Data  for  Ottawa  sand,  crushed 
quartz  sand,  and  a  graded  sand  are  analyzed  to  develop  a  model  for  the 
effect  of  initial  void  ratio  on  low-stress  dilation  behavior  of  cohe¬ 
sionless  materials;  and  critical  state  behavior  is  discussed.  Test 
results  showing  the  effect  of  fabric  anisotropy  on  dilation  behavior  are 
presented  and  a  test  on  sand-cement  indicates  the  role  of  contact  cohe¬ 
sion  in  the  physics  of  soil  strength. 

Measurement  of  Particle  Concentration  (Hardin,  1988c;  in  review). 
Particle  concentration  is  a  physical  property  of  particulate  materials 
that  has  great  influence  on  their  mechanical  behavior.  Void  ratio  or 
porosity  are  the  parameters  most  often  used  to  quantify  particle 
concentration.  The  determination  of  void  ratio  or  porosity  of  a  soil 
sample  is  part  of  the  procedure  for  many  ASTM  standards. 

Soil  samples  are  often  formed  by  placement  in  a  mold  with  rigid 
boundaries,  followed  by  striking  off  a  plane  surface  with  a  rigid 
straightedge.  Such  rigid  boundaries  and  plane  surfaces  interrupt  the 
packing  of  particles.  When  the  volume  of  the  sample  is  computed  from 
the  volume  of  the  rigid  mold  and  plane  surfaces,  additional  void  space 
at  the  boundaries  that  is  not  representative  of  the  void  ratio  or  poro¬ 
sity  of  a  repeating  element  of  the  packing  is  included  in  the  volume. 
The  computed  void  ratio  or  porosity  are  too  large  and  the  error  increa¬ 
ses  with  decreasing  sample  size. 

This  paper  describes  the  results  of  theoretical  and  experimental 
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studies  of  the  effect  of  rigid  boundaries  on  void  ratio  or  porosity 
measurement.  An  equation  is  presented  that  relates  the  true  void  ratio 
to  measured  void  ratio,  accounting  for  mold  size  and  shape,  and  particle  fl 

size  distribution.  The  relationship  is  shown  to  be  nearly  independent 
of  particle  shape  and  true  density. 

Elasticity  of  Particulate  Materials  (Hardin  and  Blandford,  1988a;  in  | 

review).  -  A  major  difficulty  in  formulating  elastic  constitutive  equa¬ 
tions  for  materials  composed  of  discrete  solid  particles  is  that  the 
stress-strain  behavior  is  rarely  exclusively  elastic.  Deformation  of  a 
particulate  material  almost  always  involves  slippage  at  contacts  between  | 

particles  as  well  as  elastic  deformation  of  individual  particles.  Elas¬ 
tic  deformation  is  relatively  small  and  is  often  obscured  by  plastic 
deformation  resulting  from  slippage,  rearrangement  and  crushing  of 
particles.  The  problem  is  to  isolate,  measure  and  model  the  elastic  I 

behavior  of  particulate  materials  in  the  presence  of  pervasive  yielding. 

The  elastic  behavior  of  classical  elasto-plastic  materials  is  iso¬ 
lated  by  unloading;  stress  paths  inside  the  yield  surface  are  assumed  to 
produce  elastic  strains.  But  particulate  materials  yield  during  un-  I 

loading;  consequently,  exclusively  elastic  behavior  is  restricted  to 
infinitesimal  increments  of  unloading.  Therefore,  differential  elastic 
constitutive  equations  for  particulate  materials  are  formulated  using 
measurements  of  elastic  stiffness  for  specimens  subjected  to  a  variety  ( 

of  ambient  states  of  stress. 

Previously,  Hardin  (1978)  developed  three-dimensional  elastic  dif¬ 
ferential  constitutive  equations  for  soils.  However,  these  equations 
(including  modifications  by  Hardin  1980)  were  not  formulated  to  assure  | 

objectivity.  The  primary  pjurpose  of  the  reformulation  presented  in  the 
paper  is  to  satisfy  objectivity.  Another  purpose  is  to  incorporate  the 
discovery  by  Roesler  (1979)  that  elastic  shear  stiffness  is  independent 
of  the  stress  normal  to  the  plane  of  shear.  | 

There  are  at  least  two  important  uses  for  the  proposed  elastic 
constitutive  equations.  They  may  be  used  to  define  the  elastic  strain 
increment  for  incremental  elasto-plastic  analysis;  accurate  definition 
of  elastic  strains  is  important  since  they  are  necessary  for  inverting  | 

the  elasto-plastic  stress-strain  relations.  Secondly,  the  equations 
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provide  a  three-dimensional  framework  for  defining  the  effects  of  state 
of  stress  and  soil  density  on  wave  propagation  velocities.  The  equa¬ 
tions  for  various  elastic  moduli  (shear  modulus,  constrained  modulus, 
etc.)  are  all  defined  by  the  three-dimensional  equations. 

Itepers  in  Preparation  -  Three  papers  are  currently  being  pursued  for 
possible  publication.  The  three  papers  are  tentatively  titled:  "Tri- 
axial  Compression  for  Cohesionless  Soils"  by  Hardin  (1988d),  "3D 
Constitutive  Equations  for  Cohesionless  Soils"  by  Hardin  and  Blandford 
(1988b),  and  "  3D  Elasto-Plastic  Finite  Element  Analysis  for  Cohesion¬ 
less  Soils"  by  Blandford  and  Hardin  (1988).  A  brief  description  of  each 
of  these  three  papers  is  presented  in  the  paragraphs  below.  Upon 
completion  of  the  papers,  copies  will  be  sent  to  the  APOSR. 

Triaxial  compression  for  cohesionless  soils  (Hardin,  1988d) 
presents^the  formulation  and  verification  of  a  model  for  the  stress- 
strain  behavior  of  cohesionless  soils  in  triaxial  compression.  The 
model  includes  work  softening  and  accurately  represents  the  volumetric 
strain  behavior.  Results  of  hundreds  of  drains'  triaxial  compression 
tests  with  volume  change  measurements  available  in  the  literature  have 
been  incorporated  into  a  computerized  database.  These  test  results  have 
been  used  to  formulate  the  model.  Hie  developed  3D  constitutive  equa¬ 
tions  use  the  triaxial  compression  model  to  define  the  Class  1  hardening 
function . 

Use  of  a  triaxial  compression  model  to  simulate  results  for  six 
tests  on  a  crushed  cal  cite  sand  is  illustrated  in  Fig.  4.  The  six 
samples  range  from  loose  to  dense  with  initial  relative  densities  D^  = 
0.21  to  0.72,  initial  void  ratios  e.  =  0.661  to  1.027  and  ol/p  =  2.04 

X  w  A 

for  all  tests.  Hie  measured  variations  of  principal  stress  ratio  and 
volumetric  strain  with  axial  strain  are  shown  in  Figs.  4a  and  4b  and  the 
corresponding  curves  computed  from  the  model  are  given  in  Figs.  4c  and 
4d. 

The  paper  ”3D  Constitutive  Equations  for  Cohesionless  Soils"  by 
Hardin  and  Blandford  (1988b)  presents  the  mathematical  formulation  of  a 
complete  model  for  monotonic  loading  of  cohesionless  soils.  Hie  objec¬ 
tive  is  to  present  a  complete  constitutive  equation  algorithm  that  can 
be  incorporated  into  finite  element  analysis  programs  with  minimimi 
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effort.  The  paper  includes  many  of  the  equations  presented  in  the  "Soil 
Elasticity"  and  "Soil  Plasticity"  sections  of  this  report.  However,  the 
equations  axe  arranged  in  the  order  and  format  necessary  for  computa¬ 
tion. 

Blandford  and  Hardin  (1988)  present  a  three-dimensional  finite 
element  implementation  for  the  constitutive  equations  of  Hardin  and 
Blandford  (1988a,  b) .  Included  in  this  effort  is  an  initial  stress 
solution  algorithm  for  solving  the  nonlinear  algebraic  equations  based 
on  a  second-order  predictor-corrector  scheme.  Computational  details 
associated  with  calculating  the  void  ratio,  load  step  scaling  at  the 
break  point  between  pre-peak  and  post-peak  (work  softening  behavior) 
Class  1  behavior  and  the  variable  load  path  algorithm  are  presented. 
Sample  experimental  problems  are  simulated  demonstrating  the  correctness 
of  the  finite  element  implementation  and  the  results  for  a  plane  strain 
footing  problem  (boundary  value  problem)  are  also  presented. 

THREE-DIMENSIONAL  ELASTO-PLASTIC  OOHESICJNLESS  SOIL  ANALYSIS 

This  portion  of  the  report  presents  the  current  elasto-plastic 
constitutive  equations  and  their  implementation  into  a  three-dimensional 
finite  element  code  EPSAP  ( Elasto-Plastic  Soil  Analysis  Program).  The 
elasto-plastic  constitutive  equations  are  being  investigated  for 
publication  (Hardin,  1988d;  Hardin  and  Blandford,  1988b  and  Blandford 
and  Hardin,  1988).  However,  the  equations  have  been  sufficiently 
developed  to  pursue  their  implementation  into  EPSAP  as  well  as  perform 
some  numerical  experiments. 

EPSAP  has  been  developed  to  solve  three-dimensional  cohesionless 
soil  problems  based  on  natural  (true)  stress-strain  elastic  and  plastic 
constitutive  equations  presented  in  Hardin  and  Blandford  ( 1988a, b). 
EPSAP  includes  a  linear  line  element  (Fig.  5a)  to  model  line  traction 
loads,  a  four  node  quadrilateral  surface  traction  load  element  (Fig. 
5b),  the  modified  eight  node  hexahedron  volume  element  of  Wilson  et  al . 
(1973)  including  the  patch  test  modifications  of  Taylor  et  al.  (1976) 
(Fig.  5c)  and  an  infinite  element  from  Marques  and  Owen  (1984)  which  is 
based  on  the  mapping  technique  of  Zienkiewicz  et  al.  (1981,  1983)  (Fig. 
5d) .  Solution  of  the  incremental  elasto-plastic  finite  element 
equations  is  based  on  an  initial  stress  formulation  which  utilizes  a 
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(a)  Two  Node  Line 
Load  Element 


V 

(b)  Four  Node  Quadrilateral 
Surface  Load  Element 


\ 
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c 

(c)  Modified  Eight  Node 
Hexahedron  Element 
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I 

I 

I 

i 

c 

(d)  Isoparametric  Eight  Node 
Infinite  Element 


Fig.  5.  Finite  Element  Library 


predictor-midpoint  iterative  corrector  scheme  to  satisfy  equilibrium.  A 
natural  (true)  stress  and  strain  finite  element  representation  is 
obtained  by  updating  the  nodal  coordinates  during  the  analysis  (Hsu, 
1986).  These  items  are  discussed  in  the  following  sections  along  with 
some  of  the  details  for  implementing  the  elastic  and  plastic  constitu¬ 
tive  models  and  the  presentation  of  some  standard  nunerical  procedure 
verification  problems. 


SOIL  ELASTICITY 

Hardin  and  Blandford  (1988a)  concluded  that  elasticity  modeling  for 
particulate  materials  needs  to  include:  (1)  geometry  of  the  assemblage 
(fabric);  (2)  elastic  properties  of  individual  particles;  (3)  stress 
history;  and  (4)  current  state  of  stress.  These  effects  were  included 
in  the  developed  constitutive  equations  by  a  scaler  function  and  three 
matrices  representing  different  parts  of  the  constitutive  behavior  for 
cohesionless  soils,  i.e. 


(dee)  =  ^  (S  ]  [Sf]  [S  ]  (do'} 
i— n  v  l  o 

Pa 

=  (S)  (do’}  (1) 

2 

where  F(e)  =  0.3  +  0.7e  approximates  the  effect  of  void  ratio  e  which 
changes  with  subsequent  stress  history;  p  is  atmospheric  pressure;  n  is 

&  fjy 

6  6606  6  6  l 

a  power  stress  coefficient;  (de  }  =  Lde  de  de  dy  dy  dy  J  is  the 
^  x  y  z  xy  yz  zx  T 

incremental  elastic  strain  vector;  (do')  =  Ido'  do’  do'  dT  dT  dT  j 

x  y  z  xy  yz  zx 

is  the  effective  incremental  stress  vector;  [S^]  is  isotropic 

elastic  Poisson  ratio  matrix  (Table  3);  (S^.]  is  the  reference  soil 

fabric  matrix  (Table  4)  which  is  defined  to  account  for  fabric  geometry, 

elastic  particle  properties,  and  stress  history  prior  to  measurement  of 

reference  fabric;  [S  ]  =  [R  ]  [Z]  [R  ]  is  the  transformed  stress- 

o  o  o 

compliance  matrix  (Table  5)  which  is  used  to  model  the  effect  of  state 

of  stress  on  the  constitutive  coefficients  including  changes  in  soil 
fabric  that  result  from  stress  history  following  the  measurement  of 
reference  fabric  (i.e.  stress  history  effects  not  accounted  for  by 
F(e));  [S]  is  the  nonlinear  elastic  material  compliance  matrix;  (  ) 
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Table  3.  Isotropic  Elastic  Poisson  Ratio  Matrix  [S^] 
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Table  4.  Reference  Soil  Fabric  Matrix  [S„] 

(x,y,z  -  Principal  Fabric  Directions) 
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Table  5.  Stress  Compliance  Matrices  (  [S^]  =  [R^]  *  [E]  [R^]  ) 
(a1  1  a 2  1  o3;  Principal  Stresses) 
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Table  5 .  Continued 
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*1»  tru  ,  n^  -  principal  stress  direction  cosines 
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symbolize  column  vector;  L  J  represent  row  vector;  [  ]  iB  used  to 

designate  a  matrix;  and  superscript  T  denotes  matrix  transpose.  Stress 
invariants  (Table  6)  have  been  used  in  formulating  the  stress-compliance 
matrix,  in  order  to  make  the  resulting  strain-stress  relation 
independent  of  coordinate  transformation  (objective). 

SOIL  PLASTICITY 

There  are  important  differences  between  classical  plasticity  and 
soil  plasticity.  In  classical  plasticity  the  direction  of  the  plastic 
strain  increment  vector  is  assumed  to  depend  on  the  state  of  stress  but 
is  considered  independent  of  the  stress  increment  vector  direction. 
This  is  not  the  case  for  soils,  as  can  be  illustrated  by  considering  the 
strain  increments  resulting  from  triaxial  compression  (TC)  and  isotropic 
compression  (IC)  stress  increments.  Figure  6(a)  shows  these  two  princi¬ 
pal  stress  increments  added  to  an  initial  isotropic  state  of  stress. 
Figure  6(b)  shows  the  resulting  principal  strain  increments.  The  two 
strain  increment  vectors  do  not  coincide.  Hence  the  direction  of  the 
strain  increment  vector  is  not  uniquely  determined  by  the  state  of  ef¬ 
fective  stress.  This  peculiar  behavior  of  soils  arises  due  to  the 
particulate  composition  of  soil  materials.  Soil  stress-strain  behavior 
is  controlled  to  a  great  extent  by  the  kinematics  of  particle  movement. 

Hardin  (1978)  shows  that  the  elements  required  to  define  soil 
plasticity  are 

*  a  definition  of  loading, 

*  plastic  potential  functions  for  the  different  classes  of 
loading,  and 

*  hardening  functions  for  the  classes  of  loading. 

The  definition  of  loading  is  used  to  specify  different  classes  of  load¬ 
ing  for  different  stress  increment  vector  directions,  and  to  specify 
whether  a  given  stress  increment  constitutes  loading  or  unloading. 
Hence,  the  definition  of  loading  and  plastic  potential  functions  define 
those  aspects  of  the  model  that  are  usually  defined  by  the  yield  func¬ 
tion  in  classical  plasticity. 

Hardin  (1978)  concluded  that  the  behavior  of  soils  can  be  repre¬ 
sented,  with  sufficient  accuracy  for  most  purposes,  by  considering  only 
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two  fundamental  classes  of  loading.  These  two  classes  are  called  Class 

1  and  Class  2  loadings.  The  deviatoric  component  of  stress  has  the 
major  influence  on  the  deformation  resulting  from  Class  1  loading, 
whereas  the  isotropic  component  of  stress  is  more  significant  for  Class 

2  loading.  A  stress  increment  may  be  exclusively  Class  1  or  Class  2,  or 
it  may  involve  a  combination  of  the  two.  Separation  of  the  two  classes 
of  behavior  is  based  on  particle  kinematics  (not  strictly  deviatoric 
versus  isotropic  components ) .  Class  1  stress  increments  produce  the 
type  of  particle  kinematics  to  which  stress-di latency  is  applicable. 
Stress-di latency  theory  is  used  to  define  the  Class  1  plastic  potential 
function.  Stress-di latency  theory  is  not  directly  applicable  to  Class  2 
loadings.  Consequently,  two  plastic  potential  functions  and  two  har¬ 
dening  functions,  corresponding  to  the  two  classes  of  loading,  are 


needed  to  define  the  plastic  strains  in  soils. 

The  plastic  strain  increments,  d€?.,  for  each  class  of  loading  can 
be  represented  as 


(de?.).  =  d>£ 
ij  k  k 


h  H 

where  dX^  and  dX^  are  the  incremental  hardening  functions  for  Class  1 
and  Class  2  loadings,  respectively;  gj,  g ^  represent  the  Class  1  and 
Class  2  plastic  potential  functions  (e.g.  Fig.  7);  and  of  .  are  the 
effective  stress  components.  The  total  plastic  strain  increment  is 
defined  as 


dcPj  =<!-*>  *  *  (d£f\)2  ,3) 

The  values  of  4<  (which  represents  the  level  of  class  participation)  are 
given  in  Fig.  8  for  different  stress  increment  paths  from  the  state  of 
effective  stress  defined  by  point  A.  Figure  8  constitutes  the  loading 
definition. 

A  simplified  explanation  of  the  difference  between  the  behaviors 
for  Class  1  and  Class  2  stress  increments  is  that  Class  1  stress 
increments  (e.g.  triaxial  compression,  plane  strain  compression  and 
simple  shear)  approach  the  strength  envelope  (failure);  whereas,  Class  2 
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(a)  Stress  Increments  (b)  Plastic  Strain  Increments 


Fig.  6.  Comparison  of  Plastic  Strain  Increments  for  TC  (Triaxial  Compression) 
and  IC  (Isotropic  Compression;  ISL  -  Isotropic  Stress  Line) 


loading  stress  increments  (e.g.  one-dimensional  strain  and  isotropic 
compression)  are  concave  toward  the  stress  axis  and  exhibit  stiffer 
behavior  with  increasing  load.  Ibis  is  reflected  in  the  different 
shapes  of  the  hardening  functions  for  Class  1  and  Class  2  loadings  as 
shown  in  Fig.  9. 


Strength  Model  -  An  understanding  of  the  strength  of  soils  is  essential 
to  the  formulation  of  general-purpose  constitutive  models  for  soils. 
Failure  is  a  limiting  case  and  should  be  defined  first,  so  that  the  form 
of  equations  adopted  for  various  elements  of  the  constitutive  model  (the 
Class  1  hardening  function  in  particular)  will  behave  properly  as  limit¬ 
ing  conditions  are  approached.  The  Mohr-Coulomb  failure  theory  is  used 
by  most  geotechical  engineers  to  model  the  strength  of  soils.  However, 
a  more  fundamental  approach  that  considers  the  particulate  nature  of 
soils  has  been  developed  by  Hardin  (1985b)  and  will  be  used  to  represent 
cohesionless  soil  strength. 

There  are  two  fundamental  mechanisms  that  govern  the  deformation  of 

particulate  materials:  (1)  bonding  at  particle  contacts;  and  (2)  the 

kinematics  of  particle  movement  within  an  element  of  deforming  soil.  In 

formulating  constitutive  equations  for  soils,  the  two  parameters  that 

represent  the  two  basic  particulate  mechanisms  are:  (1)  the  bonding 

obliquity  angle  <t>  ;  and  (2)  the  maximum  rate  of  dilation  d  .  Soil 
o  max 

strength  in  terms  of  effective  stress  has  been  defined  by  Hardin  (1985b) 
using  these  two  parameters.  The  strength  for  cohesionless  soils  in 
triaxial  compression  is  defined  by 


Triaxial  Compression  (TC) 
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Fig.  8.  Definition  of  Loading 


3-D  STRAIN  PARAMETER  3-D  STRAIN  PARAMETER 


Fig.  9.  Class  1  and  Class  2  Hardening  Function  Schematics 
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(4c) 


sin  <J>  =  k_  (  =•  -  <fr  )  tan  <t> 

CV  I  C  p  JJ 

+  (  1  -  k„)  sin  <t> 
i  p 


K  .  = 


1  sin  ♦ 


min  1  -  sin  4> 


(4d) 


1  +  sin  4> 


cv 


cvTC  1  -  sin  4> 


cv 


RmaxTC  =  RcvTC  +  (2  Kmin  ~  ^vTC*  dmaxTC 


(4e) 

(4f) 


where  d  is  the  value  of  d  at  a l  =  0;  d  is  a  parameter  related  to 

oo  maxTC  3  n 

the  maximum  negative  rate  of  dilation;  o’  is  the  dilation  reference 

a 


stress;  <?£  is  the  reference  stress  that  defines  the  rate  at  which  <t> 


decreases  with  increasing  ol  (*  10  p  );  4>  is  the  mineral  friction  angle 

j  a  p 

which  equals  the  bonding  obliquity  angle  for  cohesionless  soil;  <t>  is 

pOo 

the  value  of  4>  for  oi  =  0;  <t>  is  the  critical  state  angle  of  shearing 
p  J  cv 

resistance;  r  is  a  parameter  which  defines  the  effect  of  confinement  on 
o 

4>  ;  k_  is  the  critical  state  strength  coefficient  (-  0.60);  K  .  is  the 
p  t  min 

minimum  energy  transmission  ratio;  RcvTC  is  the  critical  state  principal 
stress  ratio  for  triaxial  compression  loading;  and  Rnax^  is  the  maximun 
principal  stress  ratio  for  triaxial  compression  loading. 

For  an  arbritrary  computation  stress  path,  the  strength  equations 
include  (Hardin,  1985b): 


Computation  Path  (CP) 


m' 


b  =  i  (1  -  X3  tan  e) 


1  ♦ 


(5a) 

(5b) 


F(b)  =  4b®  (1  -  b) 


R  =  R  ♦  2<R  -  K  .  )  F(b)  d 

max  maxTC  cvTC  mm  maxTC 


(5c) 
( 5d) 


where  b  is  a  parameter  which  defines  the  intermediate  principal  stress; 
m’  is  a  power  coefficient  for  the  interpolation  function  F(b);  0  is  the 
polar  principal  stress  invariant  (see  Table  6);  and  is  the  maximum 
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principal  stress  ratio.  Equations  5  have  generalized  Eqs.  4  with 
respect  to  stress  path  through  the  interpolation  function  F(b).  The 
interpolation  function  accounts  for  the  effect  of  stress  path  on  energy 
transmission  and  dissipation  at  the  peak  principal  stress  ratio.  F(b) 
gives  identical  strengths  in  triaxial  compression  (b=0)  and  extension 
(b=l).  It  has  also  been  chosen  so  that  the  value  of  b  =  (o^  -  - 
Og)  for  plane  strain  matches  test  results . 


Class  1  Hardening  and  Softening  -  A  Class  1  hardening/softening  func¬ 
tion  defines  the  relationship  between  two  scaler  quantities  representing 
3D  states  of  stress  and  plastic  strain,  respectively.  The  principal 
stress  ratio  R  (=o|/og)  has  been  used  by  Hardin  and  Blandford  (1988b)  to 
represent  state  of  stress  and  Class  1  plastic  work  represents  the 
state  of  plastic  strain.  Development  of  the  hardening/softening 
function  for  general  Class  1  loading  requires  that  a  relationship 
between  R  and  that  is  independent  of  Class  1  stress  paths  be 
formulated. 

Relationships  between  the  normalized  triaxial  compression  stress 
ratio  r^,  and  the  normalized  triaxial  compression  plastic  work  are 


pre-peak  behavior: 


$,P  -  1  -  ( 1  -  ra  )  ^ 
TC  u  TC' 


post-peak  behavior: 


^TC  =  1  +  ®TC 
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-  1 


TC  R  __  -  1 
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RcvTC  “  1 
W  "  RmaxTC  '  1 


(6a) 

(6b) 

(6c) 

(6d) 


in  which  a,  0,  a^,  and  m^,  are  pre-peak  and  post-peak  shape  parameters. 
The  normalized  stress  ratio  and  plastic  work  represent  a  relationship 
which  is  invariant  for  Class  1  stress  paths  (Hardin  and  Blandford, 
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1988b) . 

The  normalized  stress  ratio,  r,  and  normalized  plastic  work,  W^, 
for  an  arbitrary  Class  1  stress  path  are  defined  as  (Hardin  and 
Blandford,  1988b) 


r  = 


R  -  1 

R  -  1 
max 


(7a) 


W?=  f. 


3  J 


where 


'1  W, 


(7b) 


R  =  R  __  +  [R  -R  „]  [1  -  F(r)  ] 
max  maxTC  max  maxTC 


F  ( r )  =  4rn  (1  -  r11 ) 
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log  (0.5) 


log 


RcvTC  ~  1  1 
RmaxTC  “  1 


in  which  f^  is  a  constant  and  depends  on  d^^^  and  the  initial  void 
ratio  e  .  Equation  7a  is  iteratively  evaluated  until  r  converges  and 


R 


max 


is  initially  set  equal  to 


R  calculated 
max 


in  the  strength 


alogorithm  (Eqs.  4  and  5). 

Since  the  relationship  between  r  and  VfP  is  unique  (i.e.  independent 
of  the  Class  1  stress  path),  the  following  relationships  hold  (Hardin 
and  Blandford,  1988b) 


rTC  =  r  (8a) 

wj  =  WP,  (8b) 

Equations  6-8  are  used  to  express  invariant  plastic  stiffness  rela¬ 
tionships  in  terms  of  triaxial  compression  relationships.  Baaed  on  the 
triaxial  compression  developments,  EPSAP  simply  calculates  r  and  wP  = 
values  to  be  used  for  calculating  the  Class  1  plastic  work  in  Eq. 
7b.  The  plastic  work  of  Eq.  7b  is  used  to  calcuate  the  incremental 
plastic  work  which  is  used  to  calculate  plastic  strains  for  arbritrary 
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Class  1  stress  paths  as  discussed  following  the  next  section. 

Class  2  Hardening  -  The  Class  2  hardening  function  also  defines  the 
relationship  between  two  scaler  quantities  representing  3D  states  of 
stress  and  plastic  strain,  respectively.  Normalized  mean  principal 
stress  o'  is  used  to  represent  the  state  of  stress  and  the  Class  2 

ID 

normalized  plastic  work  represents  the  state  of  plastic  strain.  A 
preliminary  relationship  between  o^  and  that  is  independent  of  stress 
path  for  Class  2  loadings  in  terms  of  one-dimensional  strain  (see  Tables 
7,  8  and  9  for  the  one-dimensional  strain  model)  is  (Hardin  and 
Blandford,  1988b  are  developing  an  analytic  represenation  which  will  be 
similar  to  the  Class  1  stress  path  invariant  relationship) 


^m^lD  =  ^m^CP 


(W2 *CP 


W! 


ID 


(9a) 

(9b) 


Normalized  Class  2  plastic  work  is  defined  as  follows 


(R 


max 


(R 


max 


(10) 


where  is  a  constant  and  o’  s  o’/p  .  Based  on  the  developed  one- 
2  mm  a 

dimensional  strain  relationships,  EPSAP  simply  reads  in  a  table  of  o^ 
and  values  which  are  calculated  based  on  the  algorithms  given  in 
Tables  7-9.  Cubic  spline  interpolation  based  on  the  is  used  to 
calculate  the  normalize  class  2  plastic  work.  The  class  2  plastic 
work,  Wg ,  is  calculated  from  Eq.  10  which  is  used  to  calculate  the 
increment  in  plastic  work.  Calculation  of  the  Class  1  and  Class  2 
plastic  strain  increments  depends  on  the  incremental  Class  1  and  Class  2 
plastic  work  and  the  gradients  of  the  plastic  potential  surface. 
Details  of  these  calculations  are  provided  in  the  following  section  for 
arbritrary  stress  paths. 


Class  1  and  Class  2  Plastic  Strain  Computations  -  Equation  2  shows  that 


Table  7.  Algorithm  for  lD-Strain  Total  Work  Potential 
(notation  is  defined  in  Table  9) 
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Note:  The  constrained  modulus  M  =  - 


Table  8.  Algorithm  for  lD-Strain  Elastic  Work  Potential  and  Plastic 
Principal  Strain  Increment  Ratios  (notation  is  defined  in 
Table  9) 


Given:  v,  n,  Sj,  S2,  Sg,  KqI  p&,  (o{)1D»  a nd  e1D  and  w1D<Table  7) 

where  the  lD-strain  test  principal  stress  axes  coincide  with 

reference  fabric  principal  axes  for  the  soil,  [S  ]  =  [£], 

S  =  S, ,  S  =  Sof  S  =  S_;  and  K  1  1  ° 

x  1  y  2  z  3  o 


0 

Compute:  =  lD-strain  elastic  work  potential  at  (opipi  and 


(dC?/d€^)^D,  i  =  2,  3,  lD-strain  plastic  principal 
1  strain  increment  ratios 


Step 


!.  cj  *  ^  -  v  K 
1-n 


1-n 

o 


I  ♦  i 

S2  S3J 


K 

2*  ‘  ^ 
K1_n 

3-  c5 =  v 


.  K1_n  ' 
1  +  — 

S1  S3 
,  K1_n 

I  +  -O— 
S1  S2 


4.  c*  =  c^  +  Kq  (c£  +  c^) 


5.  kfo  =  (0.3  +  0.7  e^D) 


(aPlD 


,1-n 


6.  «1D  =  c*  k|D 


rd€e  1 
_ x 

de. 


cl 


t  k!o 


w 


ID 


ID 


i  =  1,  2,  3 


8. 


*? 

*? 


dil 

de- 


ID 


ID 


1  - 


de, 


i  =  2,  3 
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Table  9.  Algorithm  for  lD-Strain  Normalized  Stress  and  Plastic  Work 


0 

Given:  (oj)1£),  Kq,  X2»  «1D  and  as  functions  of  oj  fro®  Tables  7 
and  8;  and  p 

8L 

Compute:  =  normalized  stress;  and  Wp^  =  normalized  plastic  work 

Step 


("PlD 

°1D  *  (^1D  =  iV2  (1  *  2  Ko> 

a 

2.  Compute  for  =  Ko(oPjjj  using  the  Eq.  4f 


'•  =  I 


^PlD 
0 


(wlD  '  W1D) 


4.  Wp 


[  -  ko  ]  r* 


,P  -X, 
ID 


1D  iRm*K  '  U 


I  Notation  for  Tables  7-9 

(o^iD  =  major  principal  stress  for  lD-strain  loading; 
v  =  elastic  Poisson  ratio; 

Si  =  principal  direction  elastic  stiffness  coefficients  (isl,2,3); 

|  n  =  stress  power  coefficient; 

p  =  void  ratio/total  work  potential  power  coefficient; 

S„„  =  maximum  dimensionless  stiffness  coefficient  for  lD-strain; 

lDmax 

S.„  .  =  minimum  dimensionless  stiffness  coefficient  for  lD-strain; 

1  Drain 

a’  =  crushing  mechanism  break  point  stress; 

I  b 

o’  =  crushing  mechanism  reference  stress; 

UK 

=  power  coefficient  used  to  calculate  the  dimensionless  1D- 
strain  stiffness; 

S„_  =  dimensionless  lD-strain  stiffness; 

i  ID 

eQ  =  zero  principal  stress  void  ratio; 

e^  =  lD-strain  void  ratio; 

K  =  lateral  effective  stress  coefficient  for  ID- strain;  and 
o 

=  Class  2  plastic  work  power  coefficient. 

I  & 


1 
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the  Class  1  and  Class  2  plastic  strain  increments  can  be  evaluated  in 
terms  of  the  derivatives  of  the  plastic  potential  functions  with  respect 
I  to  the  effective  stresses.  However,  it  is  not  necessary  to  explicitly 

define  the  plastic  potential  functions  in  terms  of  the  effective  stress¬ 
es  and  then  differentiate.  Rather,  you  can  work  directly  in  terms  of 
the  gradients  of  the  plastic  potential  functions  and  express  the  plastic 
®  strain  increments  as  (Hardin  and  Blandford,  1988b) 

{dep)k  =  dxjj  (g’)k  (11) 

|  where 

18’ )k  =18’  By  8;  8^,  S’z  8^  J* 

in  which  (g’)k  is  the  vector  of  plastic  potential  function  gradients 

I  represented  in  Eq.  2. 

The  plastic  potential  function  gradients  are  a  function  of  the 
principal  plastic  strain  ratios  (Hardin  and  Blandford,  1988b).  Conse¬ 
quently,  the  first  step  in  computing  the  gradient  functions  is  to  calcu- 
|  late  the  normalized  principal  plastic  strain  ratios.  Class  1  normalized 

principal  plastic  strain  ratios  ((ipj  2  (deP/deP)jj  i  =  1,2,3)  are 
evaluated  as 


D  n  *2 

<*2>J  =  <  1  *1  >  b 

<*3>1  5  -  JO  ♦  3b^>  I 

where 

*  108  [  PTT  ) 

2  '  <>&» 

{  .  log  d/3) 

3  l0g 


(12a) 

(12b) 

(12c) 


i 
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\P  =  r 
PS 


(0.5)m  +  0.5  (1  -  r) 


The  parater  D  represents  the  stress-dilatancy  behavior  of  cohesionless 
soils  and  is  expressed  as  (Hardin  and  Blandford,  1988b) 


prepeak: 


D  = 


(1  +  dmaxTC)  R 


-  C 


max 


vl  R 


max 


fM'-iM 

-  max 


(13a) 


R  -  R 


postpeak : 


D  =  1  +  d 


cvTC 


maxTC 


^max  ^cvTC  J 


(13b) 


Note  that  the  maximum  rate  of  dilation  occurs  at  the  peak.  Parameter 
Cvj  is  used  to  define  a  parabolic  pre-peak  stress-dilatancy  relationship 
and  is  the  value  of  b  for  plane  plastic  strain. 

The  Class  2  normalized  principal  plastic  strain  ratios  are  (Hardin 
and  Blandford,  1988b) 


U{)2  -  1 


(14a) 


<if)2 


where 


<*i>l  + 


1  - 


dl 

P 


de 


l  J 


ID 


1 

R 


A 


;  i=2,3  (14b) 


log 


ti .  = 


de? 

dlf 


2  +  D 


♦  D 


ID 


log  (Kq) 


i=2,3 


K  is  the  lateral  effective  stress  coefficient  for  one-dimensional 
o 

strain  (ID);  and  (de?/de^)jp  are  the  one-dimensional  strain  principal 
plastic  strain  increment  ratios  as  given  in  Table  8  for  i  =  2,  3. 

The  plastic  potential  function  gradients  of  Eq.  11  are  calculated 
in  terms  of  the  principal  normalized  strain  ratios  as 


<g-ik  =  merl  («')k 


(15) 
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where 


fi’}k  =  li{  «2  83  0  0  °i 

and  [R£]  *  is  given  in  Table  5. 

With  the  plastic  strain  ratios  computed  from  Eq.  15,  the  plastic 
strain  increments  are  calculated  from  Eq.  11  where 


and  A  signfies  that  the  quantity  is  evaluated  as  the  midpoint  of  the 
increment.  The  plastic  strain  increments  computed  from  Eq.  11  are  used 
in  Eq.  3  to  calculate  the  total  plastic  strain  increment. 


FINITE  ELEMENT  PCRMULATION 

Finite  element  stiffness  analyses  are  based  on  generating  element 
force-displacement  equations  and  then  assembling  the  element  equations 
to  represent  the  soil  medium  equilibrium  equations.  Equilibrium  equa¬ 
tions  for  nonlinear  problems  are  typically  stated  in  incremental  form. 
Generation  of  the  incremental  finite  element  stiffness  equations  re¬ 
quires  proper  representation  of  the  total  incremental  strain  vector  {a€} 
and  the  incremental  effective  stress  vector  {act* }  where  a  symbolizes 
finite  increment.  The  total  incremental  strains  are  expressed  as  the 
sum  of  the  incremental  elastic  strain  vector  (a€  }  and  the  incremental 
plastic  strain  vector  (a€P) ,  i.e. 

{aO  =  {a€6J  +  (A€P)  (17) 

Substituting  Eq.  1  into  Eq.  17  and  solving  for  the  incremental  effective 
stress  vector  leads  to 

(ACT* 1  =  [S]"1  ((AC)  -  { A€P) ) 
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=  [C£]  {A€}  -  [CE]  UcH  (18) 

where  [C^  =  [S]_1  =  p*'n/F(e)  ([R^]-1  [E]"1  [RJ  [Sf]-1  [Sj'1)  is  the 
nonlinear  elastic  material  stiffness  matrix  and  the  inverse  matrices  can 
be  deduced  from  Tables  3-5;  and  {aC*3}  is  obtained  from  Eq.  3. 

Assuming  that  the  initial  state  of  the  soil  mediun  is  at  stress 
level  {o^},  the  incremental  equilibrium  condition  at  the  end  of  load 
step  m+1  is  expressed  in  terms  of  the  virtual  work  principle  for  small 
strains  as  (Washizu,  1968) 

|  m+1la’  J  S(Ue})  dfl  -  J  m+1ltJ  fi(Uu>)  dr  =0  (19) 


where  6  symbolizes  virtual  operator ;  {a£}  is  the  change  in  the  total 

strain  between  load  step  m  and  m+1  due  to  the  increment  in  the  tractions 

T 

{at}  ( { t}  =  It  t  t  J  ) ;  {u}  is  the  three-dimensional  displacement 

X  rjj  2 

vector  (=  lu  u  u  J  ) ;  Q  is  the  domain  (volune) ;  r.  is  the  traction 
x  y  z 

loaded  boundary  (surface  and  line);  and  the  pre-superscript  indicates 
the  load  step  with  m  being  the  previous  load  step  and  m+1  being  the  cur¬ 
rent  load  step.  Body  force  terms  have  been  neglected  in  Eq.  19  since 
they  are  assumed  to  be  included  in  the  initial  stress  vector.  The 
incremental  virtual  work  expression  is  obtained  from  Eq.  19  by  substi¬ 
tuting 

m+1{<7’}  =  m{o’}  +  U  O’) 

(20) 

m+1(t}  =  m{t)  +  (At) 


which  leads  to 


J  Lao'  J  S(Ue})  dn  -  J  UtJ  6((Au})  dr  =0 


Utilizing  the  usual  isoparametric  approach  to  finite  element  analy¬ 
sis  (see  Appendix  I),  the  incremental  displacements  and  strains  are 
related  to  the  element  displacements  as 

{Au}  =  [N]  Up} 
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U€}  =  [B]  Up} 


(22) 


where  [hj  is  the  matrix  of  shape  functions  describing  the  element  dis¬ 
placement  variations ;  (ip)  is  the  element  incremental  displacement 
vector;  and  [B]  is  the  element  strain-displacement  matrix  assuming  small 
strains  and  positive  compressive  strains.  Details  for  the  [N]  and  [B] 
matrices  are  given  in  Appendix  I.  Substituting  Eqs.  18  and  22  and  into 
Eq.  21  results  in 

[kj.]  Up}  =  m+1UJl)  (f)  +  UfP}  (23) 

where 

[k£]  =  |  [B]T  [C£]  [B]  dn 
^e 

(element  elastic  stiffness  matrix) 

m  =  J  in/  Itj)  dFj  ♦  f  In/  It.)  dra 

el  es 

(element  elastic  load  vector) 

UfP}  =  f  [B]T  [CU  UeP}  dn 

Jn  E 

e 

(element  incremental  plastic  load  vector) 

Additionally,  m+*UH)  is  the  incremental  load  step  multiplier;  sub¬ 
scripts  e,  1  and  s  signify  element,  line  and  surface,  respectively; 
matrices  [N^]  and  [N^ ]  are  given  in  Appendix  I;  and  the  other  symbols 
retain  their  previously  defined  interpretations.  Due  to  the  isopara¬ 
metric  implementation  of  the  elements  shown  in  Fig.  5,  the  element 
integrals  in  Eq.  23  are  evaluated  numerically  using  standard  two  point 
Gauss-Legendre  quadrature  in  each  coordinate  direction  ( see  Appendix  I ) . 

Accumulating  Eq.  23  for  all  the  elements  leads  to 

[Kg]  Up}  =  m+1UA)  iP}  +  UP15}  (24) 

The  matrices  of  Eq.  24  are  defined  as 
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IV  =  I  [Hj] 

(system  elastic  stiffness  matrix) 


{P}  =  (Pc)  +  S  tf) 

(system  elastic  load  vector) 

lAP33}  =  I  UfP) 

(system  incremental  plastic  load  vector) 

in  which  Z  signifies  surranation  over  the  elements  consistent  with  direct 
stiffness  assembly;  and  {P^}  is  the  applied  concentrated  force  vector. 
Condensation  of  the  internal  degrees  of  freedom  is  applied  to  the  modi¬ 
fied  hexahedron  element  prior  to  assembly  in  Eq.  24  (see  Appendix  I  for 
details).  Unlike  most  finite  element  stiffness  equations,  the  discre¬ 
tized  equations  of  Eq.  24  are  nonsynmetric.  This  nonsynmetry  in  the 
stiffnes'  equations  is  due  to  the  nonsynmetric  constitutive  relations  of 
the  elastic  constitutive  matrix,  [C^.] . 

The  direct  solution  of  Eq.  24  (solution  is  used  generically  here, 
the  nonlinear  solution  strategy  is  presented  in  the  next  section)  would 
be  enormous  due  to  the  large  number  of  simultaneous  equations.  Conse¬ 
quently,  the  direct  stiffness  assembly  and  solution  procedure  is  based 
on  the  unsymmetric  profile  (or  skyline)  algorithm  of  Taylor  (1977). 
Taylor’s  solution  algorithm  incorporates  the  symmetry  of  the  coefficient 
locations  with  respect  to  the  diagonal  plus  the  variable  band  structure 
of  the  nonsynmetric  finite  element  equations.  Taylor’s  algorith  has 
been  written  for  memory  resident  storage  and  solution.  Thus,  to  acco¬ 
modate  the  large  number  of  equations  which  are  encountered  in  solving 
typical  three-dimensional  problems,  Taylor’s  algorithm  has  been  modi¬ 
fied  for  out-of-core  (or  virtual  memory)  assembly  and  solution. 

NONLINEAR  SOLUTION  ALGORITHM 

The  incremental /iterative  solution  procedure  used  to  solve  the 
finite  element  discretized  equations  of  Eq.  24  for  three-dimensional 
elasto-plastic  cohesionless  soil  problems  is  discussed  in  this  section. 
Solution  of  the  nonlinear  equations  is  based  on  a  stepwise  linearization 
of  the  equations  via  an  incremental/iterative  predictor-midpoint  correc¬ 
tor  scheme  which  includes  geometric  updates  at  each  iteration.  Updating 
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the  geometry  provides  a  natural  (true)  stress/strain  representation  in 
the  context  of  the  small  strain  formulation  presented  in  the  previous 
section  (Hsu,  1966).  The  natural  stress/strain  iomulation  is  required 
to  be  consistent  with  the  elasto-plastic  constitutive  equations  deve¬ 
loped  in  Hardin  and  Blandford  ( 1988a, b).  This  solution  scheme  can  be 
classified  as  an  initial  stress  approach  (e.g.  Zienkiewicz,  1977)  for 
solving  elasto-plastic  problems  due  to  the  elastic  formulation  of  Eq. 
24.  However,  rather  than  solve  for  iterative  changes  in  the  dis¬ 
placements,  strains  and  stresses  (as  is  usual  with  the  initial  stress 
approach)  the  solution  for  the  full  incremental  values  will  be  generated 
based  on  the  midpoint  stiffness  properties.  Utilization  of  midpoint 
constitutive  properties  is  equivalent  to  using  a  second-order  iterative 
corrector  scheme  (Zienkiewicz,  1977). 

The  general  solution  strategy  in  the  load-displacement  space  can  be 
written  as 


m[Kg] ( 1  1}  Up}(l)  =  m+1(A  A)  m{P}(l-1)+  mUPP}(l_1)  (25) 

where  m  signifies  the  midpoint  load  level  (i.e.  m  is  midway  between  the 
previous  load  step  m  and  the  current  load  step  m+1 ) ;  and  the 
parenthetical  superscript  signifies  the  iteration  number.  Evaluation  of 
the  element  stiffness  matrices  and  load  vectors  used  in  Eq.  25  are  ob¬ 
tained  by  substituting 


Q  = 
e 


m  (i-l) 
e 


[B)  »  “[BJ*1'1* 


r  ,  ■  Vj'u 

el  el 


r  s  nyd-l) 
es  es 


(26) 


into  Eq.  24.  The  midpoint  geometry  for  iteration  i  is  calculated  by 

T 

updating  the  nodal  coordinate  vector  (x)  ( (x)  =  Lx  y  z  J  where  sub- 

n  n  n 

script  n  represents  node  point  n)  as 


m,  .  ( i )  ro,  , 
(x)  =  (x) 


(27) 


Equations  25-27  are  expressed  in  terms  of  the  midpoint  corrector 
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iteration  strategy.  The  predictor  step  is  obtained  for  i  =  1  and  m  =  m 
which  results  in 


“liy101  =  ”‘[ke]  ,  supi>i<°>  =  (oi 


Q' 


hk 

=  n  , 

W° 

e 

=  "V  .  , 

ij-fO) 

el  ’ 

es 

“(x)<0) 

m ,  , 

=  (x) 

(28) 


es 


Equation  28  shows  that  the  predictor  step  is  based  on  the  converged 
results  obtained  for  load  step  m. 

The  displacement  vector  solution  of  Eq.  25  is  written  as  the  sum  of 
the  incremental  elastic  displacement  vector  Up6}  and  the  incremental 
plastic  displacement  vector  UpPl ,  i.e. 

Up)U)  =  Up6}(l)  +  UpP}(l)  (29) 

where 

Upel(i)  =  (  S[KEJ<i-1)  )_1  S(P],i-1) 

(ApP)'1’  =  ["(Kjl*1-1*  “UP")11-1' 

In  order  to  evaluate  Eq.  29,  m+*US.)  must  be  known.  Startup  (i.e. 
the  first  load  increment)  of  the  predictor-midpoint  corrector  algorithm 
is  based  on  specifying  *UA.),  i.e.  the  initial  incremental  load  step. 
The  load  increment  is  adjusted  in  subsequent  load  steps  by  borrowing  a 
scheme  recommended  by  Ramm  (1981)  for  arc-length  algorithms,  i.e. 


m+1,  ..  m,  „ . 
US.)  =  US.) 


(30) 


where  mUS.)  is  the  m  load  step  increment,  mI  is  the  number  of  itera¬ 
tions  required  to  achieve  equilibrium  for  the  m^  load  step  and  1^  is 
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the  desired  lumber  of  iterations  (1^  *  3  -  5).  This  strategy  automa¬ 
tically  leads  to  smaller  load  steps  in  areas  of  severe  nonlinearity  and 
I  longer  lengths  when  the  response  is  nearly  linear.  An  upper  limit  of 

is  imposed  on  Eq.  30  to  reduce  potential  solution  drift  from  the 
equilibrium  curve. 

Accumulated  variable  quantities  through  iteration  i  for  load  step 
i  m+1  are  defined  as 


m+1{p)(i)  =  m{pl  +  uP)(i) 
m+1(pG}(1)  =mlpe)  t  Up')*11 
"+1(PP!(i)  =  ",{ppl  .  («,p)(i) 

m+1(='l(i)=  "V)  ♦  Uo’](i> 


I 

I 


(31) 


m+1  .£j(i) 

=  mU) 

+ 

<4£}(l) 

m+l{€ej (i) 

=  m{ee) 

+ 

WeV1’ 

m+l(£Pj(i) 

--  mup} 

l«p)(i) 

=  m(x) 

+ 

Up)*1* 

The  incremental  elastic  strains  and  effective  stresses  are  calculated  as 

IM6)*1*  rW11  Up6)'1’ 


Uo'l1*1  =  ICE(”o(l'U)l 


(32) 


Calculation  of  incremental  rather  than  iterative  strains  and  stresses  in 
Eq.  32  is  recommended  by  Key  et  al.  (1981)  for  materials  which  involve 
path  dependency.  Midpoint  effective  stresses  used  in  the  elastic  mater¬ 
ial  stiffness  matrix  are  calculated  as 


m(a’)(i>  =  “(a’}  +  \  Uo’}U)  (33) 

Termination  of  the  iterations  in  Eq.  25  is  based  on  an  internal 
energy  criterion  (Bathe  and  Cimento,  1980).  The  incremental  internal 
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energy  criterion  is 


UP1U)  °‘|P)(1-1> 
UpJ(I)  “{PI 


<  CT 


(34) 


-8  -4 

where  €„  is  the  energy  error  tolerance  (1  x  10  i  €„  i  1  x  10  );  and 
I6pj  1  =  lApJ  -  lApJ  which  is  the  iterative  change  in  the  dis¬ 
placement  vector. 

Evaluating  the  elasto-plastic  constitutive  equations  requires  an 
accurate  determination  of  the  element  stresses.  The  element  incremental 
strains  and  stresses  of  Eq.  32  are  evaluated  at  the  element  centroidal 
quadrature  point  which  is  the  optimum  location  for  elements  which  only 
contain  complete  linear  polynomial  representations  (Barlow,  1976).  Fur¬ 
thermore,  centroidal  stress  and  strain  computations  with  the  modified 
hexahedron  element  eliminate  the  necessity  of  recovering  the  condensed 
bubble  mode  displacements  since  their  contribution  is  zero  at  the 
element  centroid  (Hughes,  1987).  Further  details  on  stress  evaluations 
is  provided  in  Appendix  I. 


CONSTITUTIVE  EQUATION  OCMPITTATIONAL  DETAILS 

Natural  elastic  strain  modeling  requires  that  the  void  ratio  e  be 

upadated  to  reflect  changes  in  the  soil  volume.  The  initial  void  ratio, 

e.  (£  e  ),  is  defined  as 
l  o 


s 


where  V  .  is  the  initial  volume  of  voids;  V  is  the  solid  volume  which 
vi  s 

is  assumed  to  be  constant;  and  V.  =  V  .  +  V  is  the  initial  soil  volume. 

1  vi  s 

Since  the  initial  void  ratio  can  be  measured,  the  volume  of  solids  is 
calculated  as 


V 

s 


V. 

l 

1  +  e. 
i 


(36) 


Calculation  of  the  void  ratio  at  load  level  m,  iteration  i  is  obtained 
as 
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(37) 


me(i) 


V1*  - 

V 

s 


V 

s 


where  "V  ^  ^ 


is  the  centroidal  element  volume  for  iteration  i  corres¬ 


ponding  to  load  level  m.  Calculation  of  the  iterative  void  ratio  at  the 
element  centroid  is  consistent  with  the  calculation  of  the  constitutive 
properties.  The  updated  void  ratio  of  Eq.  37  is  used  in  elastic  soil 
function  F(e)  as  given  in  Eq.  1. 

If  an  anlysis  is  experiencing  pure  Class  1  loading,  EPSAP  prevents 
the  incremental  load  from  experiencing  sudden  loading  to  softening 
behavior  in  the  same  step.  That  is,  if  mt^r^  is  calculated  to  be 
greater  than  one  in  the  predictor  step  then  the  load  multiplier  is 
calculated  as 


1  - 


m+l-(l)  _  j 

m+l~(l)  m~ 
r  -  r 


(38) 


to  ensure  that  the  calculated  m+*r  will  equal  one  and  m+^aJL  is 

P 

calculated  from  Eq.  30.  The  scaling  represented  in  Eq.  38  has  been 
implemented  to  distinguish  between  loading  and  softening  Class  1 
behavior.  Class  1  behavior  after  the  scaling  of  Eq.  38  is  in  the 
softening  range. 

In  addition,  EPSAP  allows  for  the  input  of  arbitrary  stress  path 
loading.  This  is  accomplished  by  treating  the  applied  load  as  a  quasi- 
time-dependent  function.  The  piecewise  linear  varying  loads  are  input 
in  terms  of  the  total  magnitudes  for  a  given  load  interval.  For  each 
load  sequence,  different  increasing  load  magnitudes  can  be  specified  as 
well  as  different  initial  load  step  multipliers.  This  allows  the 
analysis  to  build  up  the  initial  state  of  stress  from  an  initial  zero 
state  of  stress,  for  example. 


VERIFICATION  OF  NUMERICAL  PROCEDURES 
A  sequence  of  analyses  is  presented  to  verify  EPSAP.  The  problems 
considered  are  shown  in  Fig.  10.  Figure  10(a)  represents  constant  b 
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stress  paths  to  locate  plane  strain  stress  states.  Figure  10(b)  repre¬ 
sents  constant  K  stress  paths  which  is  used  to  locate  one-dimensional 
strain.  Both  problem  sets  of  Fig.  10  represent  homogenous  stress/strain 
problems.  Hence,  only  a  single  hexahedron  element  discretization  for 
l/8t*1  of  the  synmetric  problem  is  required.  The  material  properties 
used  for  all  analyses  are: 


p  =  1.00 


n  =  0.50 


„  =  0.10 


S=S=S=S  =  S  =  S  =  1400 
x  y  z  xy  yz  zx 


d  =  0.20 
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o’  =  15 
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j’ 

CR 
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=  1.50 


p  =  0.50 


X„  =  0.71 


which  represent  typical  soil  material  properties  but  do  not  correspond 


to  any  particular  test  results  for  a  soil  specimen.  A  convergence 


-8  . 


tolerance  of  =  1x10  is  used  for  all  analyses. 

bt 


Considering  the  constant  b  stress  paths  (Fig.  10a),  the  principal 
stresses  are  increment  isotropically  to  o’  =  ol  =  ol  =  1  p  ,  followed  by 

1  Ld  O  & 


constant  b  loading  where  is  constant,  o|  is  increasing  and  o£  =  bo’  + 


(l-b)c^ 


Development  of  the  initial  state  of  stress  is  based  on  a 
reference  state  of  stress  of  zero,  i.e.  {o  )  =  {0},  and  isotropic 


compression  stress  increments  of  0.001  p  (a!  =  0.001)  for  ten  load 

6l 


steps  followed  by  isotropic  compression  stress  increments  of  0.01  p  (aS. 


=  0.01)  until  o!  =  ol  =  ol  =  1  p  .  This  is  followed  by  constant  b 

1  ^  o  A  j 

stress  paths  with  stress  increments  of  0.01  p  (  aJI  =  0.01).  Results 


for  the  constant  b  stress  path  analyses  are  shown  in  Fig.  11-13  for  b  = 
0  (triaxial  compression),  b  =  1/2  (intermediate)  and  b  =  1  (triaxial 


cr\  (increasing) 


^(constant) 


(a)  Constant  b  Stress  Paths  (0  <  b  <  I) 

a  j  (increasing ) 


°3  =  K  a  J 


(b)  Constant  K  Stress  Path  (1/4  <  K<  I) 


Fig.  10.  Load  Path  Definitions 
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extension).  Figure  11  shows  the  plot  of  the  major  principal  stress  (o|) 
versus  the  major  principal  elastic  strain,  c^,  for  the  three  values  of 
b.  Figure  11  shows  that  increasing  elastic  stiffness  in  the  direction 
of  major  principal  stress  increases  with  .  Furthermore,  Fig.  11  shows 
the  difference  in  elastic  loading  and  unloading  curves  that  results  from 
the  changing  void  ratio.  Figure  12  presents  the  major  principal  stress 
versus  the  major  principal  strain,  € ^ ,  for  the  same  three  values  of  b. 
Figure  12  shows  that  the  model  is  capable  of  representing  the  increased 
strength  associated  with  intermediate  b  stress  paths  (1  >  b  >  0), 
shifting  of  the  peak  strain  and  the  work  softening  behavior.  Figure  13 
shows  a  plot  of  the  negative  plastic  volumetric  strain  ( -e^  , )  versus 
the  major  principal  plastic  strain  (c^).  As  should  be  expected,  the 
triaxial  extension  analysis  exhibits  the  greatest  volumetric  change  with 
increasing  major  principal  plastic  strain.  Figure  13  shows  that  the 
soil  initially  experiences  volume  compression  which  is  then  followed  by 
volume  expansion. 

Considering  next  the  constant  K  stress  paths  (Fig.  10b),  the  major 

principal  stress  increments  are  specified  to  be  0.001  p  until  -  0.01 

s  i. 

p  ,  0.01  p  until  a\  =  0.10  p  and  0.10  p  until  a l  -  25  p  .  Principal 
a  *a  1  a  a  1  a  1 

stresses  =  K  Oj.  Results  for  the  one-dimensional  strain  (K  = 

1/2)  and  isotropic  compression  (K  =  1)  stress  paths  are  shown  in  Figs. 
14-16.  Figures  14,  15  and  16  show  plots  of  the  major  principal  stress 
versus  the  major  principal  elastic  strain;  major  principal  stress  versus 
the  major  principal  strain;  and  the  negative  plastic  volimietric  strain 
versus  the  major  principal  plastic  strain,  respectively.  These  figures 
show  that  the  computed  elastic  strains  are  nearly  equal  (Fig.  14) 
whereas  the  plastic  strains  are  much  smaller  for  the  isotropic 
compression  analysis  compared  with  the  one-dimensional  strain  analysis 
(Figs.  15  and  16).  The  decreased  plastic  behavior  associated  with 
isotropic  compression  should  be  anticipated  due  increased  confinement. 
Figures  14  and  15  clearly  show  the  concavity  of  Class  2  stress-strain 
curves  towards  the  stress  axis  and  exhibit  increased  stiffness  with 
increasing  load.  Figure  16  shows  that  for  isotropic  or  one-dimensional 
strain  loadings,  the  soil  exhibits  nearly  linear  volume  compression. 
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Plastic  Volumetric  Strain 
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SUMMARY 


This  report  has  presented  the  constitutive  equations  for  solving 
three-dimensional  elasto-plastic  cohesionless  soil  problems  for 
mono tonic  loading  via  the  finite  element  method.  The  soil  constitutive 
equations  have  been  carefully  constructed  to  incorporate  the  particulate 
nature  of  soils,  i.e.  kinematics  of  particle  movement,  particle  crushing 
and  particle  contact  bonding.  Constitutive  modeling  of  the  elastic 
component  is  based  on  the  Hertz  and  Mindlin  theories,  whereas  the 
plastic  component  is  formulated  using  stress-dilatancy  theory.  Since 
the  constitutive  equations  have  been  developed  utilizing  small-strain 
data  from  wave  propagation  and  vibration  tests  with  data  for  large 
strains,  the  coefficients  used  in  the  constitutive  equations  are  valid 
for  a  wide  range  of  loading  conditions. 

The  focus  of  this  research  effort  has  been  to  formulate  soil 
constitutive  equations  that  realistically  represent  the  particulate 
behavior  of  soils  This  is  being  achieved  by  using  constitutive  models 
that  have  coefficients  that  are  physically  understandable,  are  nearly 
|  constant  for  a  wide  variety  of  conditions,  and  includes  initial  efforts 

towards  understanding  the  effects  of  sample  disturbance  on  the  coeffi¬ 
cients.  Combining  the  physically  realistic  soil  constitutive  equations 
being  developed  in  this  research  with  finite  elements  that  accurately 
|  represent  boundary  conditions  and  displacement  variations  provides  a 

method  of  analysis  that  has  a  good  chance  of  predicting  field  perfor¬ 
mance. 


RECCMdKNDATIONS  FOR  Fl/WRE  RESEARCH 
Though  a  great  deal  has  been  achieved  in  understanding  and  quanti¬ 
fying  the  three-dimensional  elasto-plastic  behavior  of  soils,  much  work 
remains  to  be  completed  for  the  wide  variety  of  problems  encountered  in 
geotechnical  engineering.  Some  of  the  areas  requiring  additional  inves¬ 
tigation  include:  elasto-plastic  cohesive  soil  behavior,  cyclic  loading 
behavior  for  both  cohesionless  and  cohesive  soils,  rotation  of  principal 
stress  and  stress  increments  on  the  behavior  of  particulate  materials, 
and  inherent  and  stress  induced  plastic  anisotropy. 

Additional  finite  element  work  includes  implementation  and  testing 


59 


of  the  constitutive  equation  developments  and  the  utilization  and  veri¬ 
fication  of  EPSAP  for  a  variety  of  boundary  value  problems  encountered 
|  in  geotechnical  engineering.  This  phase  will  also  require  modifying  the 

constitutive  equations  for  the  special  cases  of  plane  strain  and  axi- 
symmetric  behavior.  Furthermore,  investigation  and  implementation  of 
large  strain  behavior  must  still  be  completed.  The  use  of  the  developed 
j  natural  strain  approach  may  not  be  computationally  feasible  for  some 

soil  static  problems. 

Another  major  activity  which  should  be  undertaken  is  the  develop¬ 
ment  of  a  computer  simulation  problem  to  investigate  the  behavior  of 
|  particle  packings.  A  preliminary  effort  in  this  regard  has  been  comple¬ 

ted  during  the  current  contract  as  outlined  in  the  next  section. 

COMPUTER  SIMULATION  OF  PARTICLE  BEHAVIOR 

When  attempting  to  explain  seme  typeB  of  soil  behavior  one  may  be 
led  to  speculate  about  the  nature  of  particle  movements  within  an  ele¬ 
ment  of  deforming  soil.  Which  contacts  within  the  element  are  experi¬ 
encing  gross  sliding  at  a  given  instant  of  deformation'1  What  is  the 
(  distribution  of  orientations  of  sliding  contacts?  Do  groups  of  parti - 

cles  slide  past  one  another,  and  if  so,  how  many  particles  are  in  the 
sliding  groupjs?  Because  it  is  difficult  to  observe  the  behavior  of 
individual  particles  and  contacts  within  a  soil  sample  in  the  labora- 
|  tory,  it  will  apparently  be  instructive  to  simulate  the  behavior  of  a 

packing  of  particles  with  a  computer  model. 

Aspects  of  soil  behavior  that  invite  questions  about  the  nature  of 
particle  movements  are:  (1)  the  effect  on  soil  deformation  of  a  change 
in  direction  of  the  stress  increment  vector  (rotation  of  stress  incre¬ 
ment  ) ;  ( 2 )  the  nature  and  development  of  inherent  and  stress -induced 
anisotropy  in  particulate  materials;  and  (3)  the  effects  of  contact 
cohesion  on  soil  deformation. 
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APPENDIX  I 


FINITE  ELEMENT  FORMULATION  DETAILS 


FINITE  ELEMENT  FORMULATION  DETAILS 
The  purpose  of  this  appendix  is  to  present  the  finite  element 
formulation  details  associated  with  the  three-dimensional  computer 
program  EPSAP  (Elasto-Plastic  Soil  Analysis  Program).  These  formulation 
details  include:  (1)  element  shape  functions  and  shape  function  matrices 
(which  includes  the  element  strain-displacement  matrices)  for  the  finite 
elements  shown  in  Fig.  5,  (2)  curvilinear  coordinate  transformations  for 
the  isoparametric  implementation,  (3)  numerical  integration  and  (4) 
stress  and  strain  calculations. 


Element  Shape  Functions  and  Shape  Function  Matrices 

Shape  functions  for  the  elements  of  Fig.  5  and  the  shape  functions 
matrices  are  explicitly  given  in  this  section.  First,  the  shape  func¬ 
tions  corresponding  the  the  line  load  element  of  Fig.  5a  are 


N*  =  (1  -  0/2 

Ng  =  (1  +0/2 


(Al) 


where  -1  <  £  i  1  is  the  nondimens ionali zed  local  coordinate  shown  in 
Fig.  5a  and  subscript  i  (i=l,2)  corresponds  to  the  two  element  node 
points.  The  matrix  of  line  element  shape  functions,  [N  ),  is  defined  as 
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Equation  A 2  shows  that  [N  ]  is  composed  of  the  two  basic  shape  functions 
(Eq.  Al ) ,  the  positioning  of  the  shape  functions  being  dictated  by  the 
nodal  x-,  y-  and  z-axis  line  traction  values.  The  line  tractions  are 
approximated  as 


i 


1 


t8- 

r  n1  o  o  i 

t*. 

X 

1 

XI 

t1 

0  N8-  0 

t1. 

7 

l 

yi 

t8- 

0  0  N8, 

t8-. 

l  z 

L  l  J 

1  Z1 

(A3) 
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where  summation  on  the  repeated  subscripts  is  implied  and  t  t  ^ ,  t^ 
are  the  nodal  line  traction  values. 

The  four  node  quadrilateral  element  (Fig.  5b)  shape  functions  are 
(e.g.  Zienkiewicz,  1977) 


(A4 ) 


in  which  -1  1  i  1  are  normalized  curvilinear  coordinates  (see  Fig. 
5b),  =  ft1*,  and  are  ±1  depending  on  node  location  i, 

e.g.  i .  -  -1  for  i=l,4  and  fl.  =  1  for  i=3,4. 

l  ’  i  ’ 

s 

The  matrix  of  surface  element  shape  functions,  [N  ],  is  defined  to 
be 


[NS]  = 


N®  0  0 

l 

0  N®  0 

l 

0  0  N® 

l 


(A5) 


in  which  the  pattern  of  Eq.  A5  is  repeated  for  i=l,2,3,4.  The  surface 
s 

tractions  (t  )  are  interpolated  in  the  same  manner  as  the  line  element 
tractions  i.e. 


r  t®  ] 

N®  0  0  1 

tS. 

X 

l 

XI 

t® 

— 

0  N®  0 

t®. 

y 

l 

yi 

t® 

0  0  N8 

t®. 

z 

l  J 

Z1  J 

(A6) 


s  s  s 

where  simulation  on  the  repeated  subscripts  is  implied  and  t^,  t^, 
are  the  nodal  surface  traction  values. 

The  element  shape  functions  for  the  basic  eight  node  hexahedron 
element  (Fig.  5c)  are  (e.g.  Zienkiewicz,  1977) 


I 


2 


NA  =  |  (1  ♦  CQ)  (1  +  ^q)  (1  +  <Q)  J  i=l,2 . 8  (A7 ) 

where  and  -1  W  i  1  ((.  =  1  for  i=2,3,6,7;  =  1  for  i  =  3,4, 
7,8;  =  1  for  i=5,6,7,8).  These  basic  eight  node  hexahedron  element 
shape  functions  are  vised  to  describe  the  element  geometry  and  the  volume 
element  shape  function  matrix,  [N].  Similar  to  the  line  and  surface 
element  shape  function  matrices,  the  volume  element  shape  function 
matrix  is  defined  as 


[N] 


N.  0  0 

l 

0  N.  0 

l 

0  ON. 

l 


(AS) 


where  i= 1,2, ... ,8. 

A  disadvantage  of  the  volume  element  shape  functions  represented  oy 
Eq.  A7  is  that  they  cannot  represent  general  linear  stress  variations. 
Wilson  et  al .  (1973)  proposed  that  the  incompatible  "bubble  mode"  func¬ 
tions: 


Nj  =  (1  -  t2) 

=  (1  -  t,2)  (A9) 

N3  =  (1  -  C2) 

which  are  zero  at  all  eight  node  points  be  added  to  element  displacement 
approximation  in  order  to  represent  linear  stress  variations.  The  addi¬ 
tion  of  the  bubble  mode  shape  functions  leads  to  the  following  element 
displacement  variations: 

u  =  N.  u.  +  Nb  a^ 

IX  J  X 

v  =  N.  v.  +  Nb  aj  (A10) 

I  1  j  y 

w  =  N.  w.  +  Nb  a^ 

II  j  z 

where  u,  v,  w  are  the  x-axis,  y-axis  and  z-axis  displacements,  respec¬ 
tively;  u^,  v^,  w^  are  the  standard  node  point  displacements;  and  a^, 
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a^,  aJ  (j=l,2,3)  are  additional  degrees  of  freedom  associated  with  the 

y  ^ 

internal  bubble  mode  shape  functions.  Utilizing  the  displacement 
approximations  of  Eq.  A10  results  in  a  modified  eight  node  hexahedron 
element . 

Equation  A10  shows  that  the  number  of  displacement  variables  for 
the  modified  eight  node  hexahedron  is  33  rather  than  the  usual  24.  Eli¬ 
mination  of  the  additional  nine  variables  (aJ,  aJ  aJ;  j=l,2,3)  is 

x  y  z 

accomplished  using  static  condensation  on  the  element  stiffness  matrix. 
Before  pursuing  static  condensation,  generation  of  the  element  strain- 
displacement  matrix  for  the  modified  eight  node  hexahedron  element  will 
be  presented. 

The  element  strain-displacement  matrix  is  obtained  by  first  de¬ 
fining  the  incremental  strain-displacement  relationships: 


d€ 

X 

i 

0 

0 

d€ 

y 

0 

_d_ 

~dy 

0 

d€ 

z 

0 

0 

d 

~dz 

du 

dy 

xy 

— 

d_ 

dy 

d_ 

dx 

0 

dv 

dw 

dy 

yz 

0 

d_ 

dz 

d 

dy 

S 

■b 

d_ 

■  dz 

0 

d_ 
dx  - 

where  the  negative  sign  for  the  direct  strain  increments  is  vised  since 
compressive  strains  are  positive.  Substituting  the  displacement  varia¬ 
tions  of  Eq.  A10  into  Eq.  All  results  in 
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-aNi/ax  o 
0  -Sti^/ay 


sti^/dy  dN./ax 
0  dNi/dz 

Sti^dz  0 


-aN./dz 


Sti^/dy 

aN./ax 


-dN^/dx 

J 


8Nb/dy 

3 


Stih/dz 

3 


-dNb/ay 

3 


aNb/ax 

3 

aNb/dz 

J 


-aNb/az 

0 


aNb/dy 

3 

aNb/ax 

J 


{de}  =  [BH]  {dp6}  +  [BB]  {da} 


(A12 ) 


where  [Bn]  is  the  strain-displacement  matrix  associated  with  the  element 

e  B 

node  point  displacements  {p  }  and  [B  ]  is  the  b train-displacement  matrix 
corresponding  the  the  bubble  mode  displacement  variables  {a} . 

The  element  elastic  stiffness  matrix  is  constructed  as 


[k-g]  =  |  [B]T  [Ce]  EB]  dn 


[bh)t  H  b 

-  1C  ]  [  [BM]  |  [B°]  ]  dfi 

[BB]T 


[BH]T  [C^J  [BH)  I  [BH]T  [O.J  [BB]  f 


EB"] ‘  [OJ  EB"}  |  [BB1T  [C£]  [BB] 


H  B 

where  [B]  =  [  [B  ]  |  [B  ]  ]  is  the  element  strain-displacement  matrix. 

Static  condensation  of  the  element  elastic  stiffness  matrix  is 
obtained  by  eliminating  the  incremental  bubble  mode  displacement 
parameters  from 


r,  HH,  HB.  ,  .  e, 
[kg  3  [k^,  ]  {dp  } 

r.  BH,  r,  BB,  .  ,  , 

[kg  1  [kg  ]  {da} 


( A14 ) 


Expanding  the  lower  partition  of  Eq.  A14  results  in 

[k^1]  {dpe}  +  [k^B]  {da}  =  {0} 

{da}  =  -  [k^V1  [k^1]  {dpe}  (A15 ) 

Substituting  Eq.  A15  into  the  upper  partition  of  Eq.  A14  leads  to 
L  [kg.  ]  -  [k^.  ]  [kg  ]  [kg  ]  ]  {dp  }  =  {df} 

[kg  ]  {dp6}  =  {df}  (A16 ) 

Equation  A16  defines  the  statically  condensed  elastic  stiffness  matrix. 

A  problem  with  the  formulation  given  by  Eq.  A13  is  that  the  modi¬ 
fied  hexahedron  element  does  not  pass  the  patch  test  unless  the  element 
geometry  is  at  least  a  parallelpiped.  Lack  of  satisfying  the  patch  test 
may  lead  to  erroneous  results.  Taylor  et  al .  (1976)  designed  a  "repair" 
for  the  two-dimensional  version  of  the  element  such  that  for  an  arbi¬ 
trary  element  geometry  the  modified  element  will  pass  the  patch  test. 
Their  repair  involves  calculating  the  bubble  mode  strain-displacement 
matrix  using  the  centroidal  (i.e.  £  =  n  =  =  0 )  evaluation  of  the 

Jacobian  matrix  and  Jacobian  determinant.  Further  details  of  this 
calculation  are  presented  in  the  curvilinear  coordinate  transformations 


I 

I 


» 


and  numerical  integration  sections. 

The  infinite  element  shown  in  Fig.  5d  is  taken  from  Marques  and 
Owen  (1984)  and  is  based  on  the  mapped  infinite  element  concept  of 
Zienkiewicz  et  al.  (1981,  1983).  Mapped  infinite  elements  are  based  on 
a  simple  mapping  technique  that  applies  to  both  modelling  of  the  geome¬ 
try  and  the  field  variable  (e.g.  displacements  in  stress  analysis). 

The  concept  underlying  mapped  infinite  elements  is  most  easily 
understood  in  the  context  of  the  one-dimensional  element  shown  in  Fig. 
A.l.  Nodes  1,  2  and  3  (node  3  is  at  infinity)  of  Fig.  A. la  are  mapped 
onto  the  parent  element  defined  by  the  local  coordinate  system  -lit! 
1  as  shown  in  Fig.  A. lb.  The  pole  (i.e.  point  of  singularity  for  the 
mapping)  position  labelled  0  in  Fig.  A. la  is  arbitrary  (though  its 
position  influences  the  accuracy  of  the  results)  and  Xq  1  x^.  Once  x^ 
is  chosen,  the  location  of  node  2  is  defined  by 

X2  =  ^X1  “  x0  < A17 ) 


Interpolation  between  the  local  and  global  coordinate  systems  is 


2 

x(£)  =  I  M .(C)  x  (A18) 

i=l 


where  the  summation  extends  over  the  finite  nodes  only  and  the  mapping 
functions  M^  are  given  by 


Mj  =  -2£ / ( 1  -  £) 

M2  =  (1  +  £ )/( 1  -  £ ) 


(A19) 


Examining  Eq.  A19  shows  that  for  £  =  -1,  0,  1  the  corresponding  x -coor¬ 
dinates  are  x^,  x^,  00 • 

Interpolation  for  the  x-axis  displacement  uses  standard  Lagrangian 
approximations  to  give 


u  =  I  Ni  u.  =  0.5£ (£-1 )  u  +  ( l-£  )  +  0.5£(£+l)  u3  (A20) 

i  =  1  1 


Solving  Eq.  A18  for  £  leads  to 


£ 


I 


1  -  2a/r 


(A21) 


I 


7 


in  which  r  denotes  the  distance  from  the  pole  to  a  general  point  within 
the  element  and  a  =  -  Xj  as  shown  in  Fig.  A. la.  Substituting  Eq.  A21 
into  Eq.  A20  results  in 

u  =  u3  +  (-Uj  +  4u2  -  3u3)  a/r  +  (2u1  -  4u2  +  2u3>  (a/r)2  (A22) 

Thus,  as  r  tends  to  infinity,  u  approaches  u3  which  is  assumed  to  be 
zero  in  the  present  implementation.  Consequently,  summation  over  the 
finite  nodes  is  all  that  is  required  in  formulating  the  finite  element 
equations.  Equation  A22  also  clearly  shows  the  role  of  the  pole  posi¬ 
tion,  0. 

The  mapping  and  displacement  shape  functions  corresponding  to  the 
eight  node  isoparametric  element  of  Fig.  5d  are  ( Marques  and  Owen,  1984) 


Mi  =  \  (i  +  CQ)  (1  +  *0>  (-f)/(l  -  C)  ;  i=l,2,3,4 
Mi  =  l  <1  +  £0)  (1  +  V  (1  +  <r)/(i  -  <r)  ;  i=5,e,7,8 


( A23 ) 


Ni  =  §  (1  +  £0)  (1  +  V  (<2  ~  °  J  i=1-2’3’4 

Ni  =  \  (1  +  £0)  (1  +  \)  (1  -  CZ)  J  i=5,6,7 ,8 


(A24 ) 


Equations  A23  and  A24  show  that  the  eight  node  isoparmetric  infinite 
element  uses  a  bilinear  approximation  in  the  finite  C  -  't1  plane  and  a 
quadratic  decay  in  the  infinite  £ -direction. 

The  infinite  element  strain-displacement  matrix  is  obtained  by 
substituting  the  shape  functions  of  Eq.  A24  into  the  incremental  strain- 
displacement  relationships  of  Eq.  All  which  results  in 
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aN./ax 

0 

0 

0 

-dN./dy 

0 

0 

0 

-aNi/az 

AN^/dy 

AN./ax 

0 

0 

a^/az 

ANi/dy 

dNi/dz 

0 

aht/ax 

(A25 ) 


The  infinite  element  elastic  stiffness  matrix  is  obtained  by 
substituting  Eq.  A25  into 


fk£]  =  J  [B]t  [ce]  [B]  dn 


CurvilinesLr  Coordinate  Transformations 

'Tie  element  matrices  and  integrals  of  the  previous  section  involve 
global  coordinates  whereas  all  the  shape  functions  have  been  expressed 
in  terms  of  local  coordinates  on  the  appropriate  parent  elements. 
Consequently,  a  transformation  of  coordinates  between  the  global  and 
local  coordinate  systems  is  required  to  formulate  the  element  stiffness 
matrices  and  load  vectors.  These  transformations  are  obtained  using 
isoparametric  finite  element  transformation  procedures  which  are  des¬ 
cribed  in  this  section. 

Basically,  an  isoparametric  finite  element  formulation  involves 
approximating  the  finite  element  geometry  (represented  by  shape  func¬ 
tions  It  )  in  the  same  manner  as  the  element  displacement  variations  in 
terms  of  the  node  point  displacements  (represented  by  shape  functions 
N..  ) ,  i.e.  set  except  for  the  infinite  element  of  Fig.  5d  in 
which  case  the  special  geometric  interpolation  functions  of  Eq.  A23  are 
used.  Details  of  the  various  transformations  required  in  EPSAP  are  pre¬ 
sented  in  the  following  paragraphs. 

The  elements  used  in  EPSAP  require  the  evaluation  of  the  arc  length 
trans  f  onnat  i  on 


dr^  =  jfu)  di 


( A26 ) 


where  J^(C>  is  the  arc  length  Jacobian,  the  surface  area  transformation 


( 
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(A27 ) 


drs  =  Jr(£»‘n)  d£ 

where  J^(£ ,fl)  is  the  surface  area  Jacobian ,  the  volume  transformation 

dfl  =  |J(£,-n,<r)|  d£  d'n  d£  (A28) 

where  |J(£»1»,£’)|  is  the  three-dimensional  Jacobian  determinant,  as  well 
as  volume  shape  function  derivatives  with  respect  to  x,  y  and  z  as  shown 
in  Eqs.  A12  and  A25. 

The  arc  length  Jacobian  of  Eq.  A26  is 

jr<o  =  1  <|r)2  +  <lr)2  +  (lf)2  <A29) 

The  finite  element  representation  of  Eq.  A29  is  obtained  by  substituting 

x  =  IM*J  (x) 

7  =  tMAJ  {y} 
z  =  IM^J  {z} 


into  Eq.  A29  which  results  in 


Jp<€ )  = 


.  [ d#  «  M d#  M 


.  .  .  ( A30) 

where  LM^-  J  =  IN^J  is  the  row  vector  of  line  element  shape  functions  (Eq. 
Al)  and  {x} ,  { y } ,  {z}  are  the  column  vectors  of  element  node  point 
coordinates.  The  sure  length  transformation  results  in  the  limits  of 
integration  being  -1  1  £  £  1. 

The  surface  area  transformation  of  Eq.  A27  is 


JpU.o) 


,dx  dy 

~  J  d£  d* 


8x  ay,  2  .dy  dz  dy  dz ,  2  . dz  dx  dz  dx,  2 

d*n  d£  '  '  d£  d-n  d*  d£  '  '  d£  d*  “  d*  d£  ' 


(A31) 


and  the  finite  element  form  is  obtained  by  substituting 
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x  =  LM8  J  {x} 
y  =  LM3  i  (y) 
z  =  LM8  J  {z> 

into  Eq.  A31;  LM8  J  =  IN8  J  is  the  row  vector  of  surface  element  shape 
functions  given  in  Eq.  A4;  and  the  other  symbols  are  as  previously 
defined. 

The  evaluation  of  the  element  volume  integrals  require  not  only  the 
volume  transformation,  but  also  requires  the  evaluation  of  x,  y  and  z 
coordinate  derivatives  as  shown  in  the  strain-displacement  matrices  of 
Eq.  A12  and  A25.  Since  the  shape  functions  are  expressed  in  terms  of 
the  nondimens ionali zed  £,  -n  and  Z  coordinates;  the  x,  y  and  z  coordinate 
derivatives  must  be  related  to  the  £ ,  Z  coordinate  system.  Using  the 
chain  rule  of  differentiation 


aN. 

i 

ox 

aN. 

_ i 

dy 

aN. 

i 

.f 

dz 

an. 

i 

d£ 

ac 

ax 

a£ 

ay 

a£ 

dz 

Ijs 

ax 

aN. 

i 

4. 

dy 

aN. 

i 

+ 

az 

aN. 

1 

a-n 

a-n 

ax 

• 

a-n 

dy 

a-n 

dz 

aN. 

i 

ax 

aN. 

i 

dy 

»i 

dz 

aN. 

i 

dZ 

dZ 

ax 
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where  N^  are  the  shape  functions  given  in  Eq.  A7  or  A24. 
A3 2  in  matrix  form  results  in 


(A32 ) 


Writing  Eq. 


!  51 

dx 

ay 
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d£ 
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dx 

ay 

a-n 

a-n 

an 
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dx 

dy 

* 

dZ 

a  z 

dz 

r  aN. 

i 

d£ 

ax 

dz 

aN. 

i 

a-n 

ay 
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dNi 
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r  3N- 
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dx 

aN. 
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(A33) 


where  IJ]  is  the  Jacobian  matrix.  The  shape  function  derivatives  with 
respect  to  x,  y  and  z  can  be  obtained  from  Eq.  A33  by  premultiplying 
both  sides  by  the  inverse  Jacobian  matrix,  [J]  ,  i.e. 
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f  aN. 

r  aN.  i 

_ 1 

dx 

d£ 

aN. 

aN. 

i 

dy 

=  [J]'1 

i 

a-n 

aN. 

aN. 

i 

i 

.  dz 

.  a^  . 

Since  the  £ ,  and  {  shape  function  derivatives  are  straightforward,  Eq. 
A33  provides  the  partial  derivatives  of  the  shape  functions  with  respect 
to  x,  y  and  z. 

The  Jacobian  matrix  of  Eq.  A33  is  more  explicitly  expressed  as 


aiMj 

d£ 

(X) 

dLMJ 

d£ 

{yl 

aLMj 

d£ 

(z) 

'  Ju 

J12 

Jl3' 

[J]  = 

dLMJ 

a-n 

(X) 

aLMj 

a-n 

{y} 

dLMJ 

a* 

(z) 

= 

J21 

J22 

J23 

dLMJ 

(X) 

dLMJ 

a<r 

(y) 

aLNy 

a<r 

(z) 

•  J31 

J32 

J33  • 

.  .  .  (A35) 


where  LMJ  =  LNJ  for  the  modified  eight  node  hexahedron  element  (see  Eq. 
A7)  and  LMJ  is  the  row  vector  of  mapping  functions  given  in  Eq.  A23  for 
the  isoparametric  eight  node  inifinte  element.  The  Jacobian  determinant 
and  inverse  Jacobian  matrix  are 


[J] 


-1 


1 

Ul 


!)  1  =  J11(J22J33 

~  J23J32)  +  J12(J 

23J31  ‘  J21J32) 

+  J13(J: 

23J31  *  J22J31) 

( A36 ) 

J22J33  "  J23J32 

J32J13  _  J12J33 

J12J23  ”  J22J13 

J31J23  _  J21J33 

J11J33  "  J13J31 

J21J13  "  J11J23 

•  J21J32  "  J31J22 

J31J12  '  J11J32 

J11J22  '  J12J21  - 

.  .  .  (A37 ) 

Substituting  Eq.  A36  into  Eq.  A28  gives  the  desired  volume  transfor¬ 
mation  for  the  hexahedron  element.  Equation  A37  gives  the  explicit  form 
of  the  inverse  Jacobian  matrix  which  is  required  in  Eq.  A34. 

As  mentioned  in  the  previous  section,  a  special  evaluation  of  the 
bubble  mode  strain-displacement  matrix  (see  Eq.  A12)  must  be  used  in 
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order  for  the  modified  eight  node  hexahedron  element  to  pass  the  patch 
test  for  an  arbitrary  element  geometry.  This  modification  involves 
evaluating  the  bubble  mode  shape  function  derivatives  as  (Taylor  et  al. , 
1976  and  Cook,  1981) 


!$ 

dx 

aNb 
_ j 

ay 

L  8z 


|J 


Ul 


IJc> 


-1 


aN~ 
_ j 

a t 

dNb 
_ J 

aNb 
_ j 


( A38 ) 


where  |J|  is  the  Jacobian  determinant  of  Eq.  A36,  | |  is  the  Jacobian 

determinant  evaluated  at  the  element  centroid  (i.e.  |J  |  =  | J(£=*>=£=0) | ) 

-1  C 
and  [J^]  is  the  inverse  Jacobian  matrix  of  Eq.  A37  evaluated  at  the 

element  centroid.  With  the  modification  of  Eq.  .*38  the  modified  eight 

node  hexahedron  element  passes  the  patch  test.  Consequently,  this 

element  converges  to  the  correct  equilibrium  solution  with  increasing 

mesh  refinement. 

With  the  curvilinear  coordinate  transformations  presented  in  this 
section,  the  isoparametric  finite  element  equations  can  be  generated. 
Unfortunately,  due  to  the  algebraic  complexity  of  the  transformations, 
closed  form  evaluation  of  the  element  integrals  is  generally  not 
possible.  Evaluation  of  the  element  integrals  is  the  topic  of  the  next 
section. 


Numerical  Integration 

Due  to  the  complications  introduced  by  isoparametric  interpolation 
into  the  element  stiffness  matrix  and  load  vector  integrals,  numerical 
integration  (quadrature)  must  be  used  to  evaluate  the  integrals.  Typi¬ 
cally,  in  finite  element  applications,  Gaussian  quadrature  formulas  are 
used  due  to  their  superior  accuracy  for  a  given  number  of  quadrature 
point  function  evaluations.  In  the  following  paragraphs;  line,  surface 
and  volume  integral  quadrature  formulas  are  presented. 

The  line  traction  integral  is  represented  as 

f  [GaU)]  JpU)  dC  (A39) 

-1  1 
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I 


where  [G^(£)]  =  [N^J^  [N®-].  Numerical  evaluation  of  Eq.  A39  is  based  on 
Gauss-Legendre  quadrature  which  is  represented  as 

I 

1  nG 

f  [G^( £ ) 3  jJ(£)  d£  J  Z  w.  [Ga(C.)]  jJ(C.)  (A40) 

-1  1  i=l  1  ill 

IG 

where  n  is  the  number  of  Gauss-Legendre  quadrature  points,  w.  is  the 

th  th  ^ 

i  weighting  coefficient  and  £ .  is  the  i  quadrature  point.  Hie  value 

G  ^ 

of  n  used  in  this  report  is  two  and  the  corresponding  Gauss-Legendre 
data  is  given  in  Table  A.l. 

I  Evaluating  the  quadrilateral  surface  load  integral  and  the  element 

volume  integrals  is  easily  obtained  by  extending  the  one-dimensional  in¬ 
tegral  of  Eq.  A40  into  the  second  and  third  dimensions,  respectively. 
Consequently,  the  surface  load  integral  for  the  quadrilateral  element 
1  becomes 


1 

-1 


[GsU,-n)3  d£  du 


n  n 

i  i 

i  =  l  j=l 


wi  fGS«i'V]  J>i’V 


(A41) 


where  [ Gs ( € , "n ) ]  =  [N5]^  [N5]  and  the  volume  integral  for  the  modified 
hexahedron  and  infinite  elements  becomes 


r1  r1  r1  v 

I  [Gv(£,n,f>]  |J(£,*»,0  I  d£  di  d<f 

J  «  J  i  J  i 


■1-1-1 


G  G  G 
n  n  n 

Z  Z  Z  w  w  w  [Gv(£ .  ,-n  )  ]  |J(C .  )  |  (A42 ) 

i  =  l  j=l  k=l  1  J  K  1  J  k  1  3  k 


in  which  [GV  (£,*>,£)]  =  [B]^  [C  ]  (B].  Hie  global  derivatives  of  the  N. 

ep  l 

shape  functions  used  in  the  element  strain-displacement  matrices  (Eqs. 
A12  and  A25)  are  calculated  using  the  transformation  of  Eq.  A34  whereas 
the  bubble  mode  shape  function  derivatives  are  obtained  using  the  trans¬ 
formation  of  Eq.  A38. 
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Stress  and  Strain  Calculations 


The  soil  elasticity  constitutive  relationships  which  are  used  to 

construct  the  elastic  constitutive  matrix,  [CL,] ,  require  an  accurate 

determination  of  the  effective  stresses  {a’}  =  lo’  o’  o'  T  7 
,p  xx  yy  zz  xy  yz 

T  J  .  Displacement  based  finite  element  analyses  typically  result  in 
zx 

accurate  displacement  approximations  but  much  less  accurate  stress 
representations.  Barlow  (1976)  investigated  the  existence  as  well  as  a 
method  of  locating  optimal  points  for  calculating  accurate  stresses  for 
displacement  based  finite  element  formulations.  His  technique  involves 
subjecting  the  finite  element  to  a  complete  polynomial  field  of  one 
order  higher  than  the  complete  polynomial  representation  included  in  the 
shape  function  approximation.  The  objective  being  to  locate  unique 
positions  within  the  element  at  which  the  stresses  have  the  same  degree 
of  accuracy  as  the  nodal  displacements.  Hence,  the  terminology  "optimum 
stress  locations." 

Barlow  (1976)  found  that  the  optimum  stress  point  is  the  element 
centroid  (i.e.  C=n=^=0)  for  the  eight  node  hexahedron.  Numerical 

experimentation  by  the  authors  has  verified  that  the  centroid  is  also  an 
optimum  location  for  the  modified  eight  node  hexahedron  element.  Con¬ 
sequently,  the  element  strains  and  stresses  used  to  generate  the  elasto- 
nlastic  constitutive  matrix  are  based  on  the  element  centroid  for  both 
the  modified  eight  node  hexahedron  element  (Fig.  5c)  and  the  eight  node 
isoparametric  infinite  element  (Fig.  5d) .  Extenting  the  centroid  stress 
evaluation  concept  to  the  infinite  element  is  based  on  the  fact  that  its 
shape  function  representation  is  a  complete  linear  polynomial  in  terms 
of  (,  m  and  f . 

Whe'1  evaluating  the  modified  hexahedron  element  stresses  at 
arbi itrary  points,  the  bubble  mode  displacement  parameters  must  be  known 
in  order  to  define  the  element  strains  (see  Eq.  A12).  Recovery  of  the 
element  bubMe  mode  parameters  are  obtained  from  Eq.  Al">  once  the  incre¬ 
mental  displacements  have  been  generated  from  the  nonlinear  solution 
procedure.  However,  since  (B  (£=n=£=0)]  is  identically  zero,  no  recov¬ 
ery  of  the  bubble  mode  displacements  is  required  (Hughes,  1987).  Thus, 
the  e'ement  incremental  strains  are  evaluated  as 

(de)  -  [BH(U“nr^=0)]  (dpe)  (A43) 

for  the  modified  eight  node  hexahedron  element  and 
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( A44 ) 


{d€ }  =  [B(£=*l=f=0) ]  {dpe} 

for  the  isoparametric  eight  node  infinite  element.  Element  incremental 
effective  stresses  are  also  evaluated  at  the  centroid  as 

{do’}  =  [Cg]  {del  (A45) 

where  [C„]  and  {d€ }  ere  both  evaluated  at  the  element  centroid, 
fc- 


Table  A. 1  Two  Point  Gauss -Legendre  Quadrature  Data 


i 

£  . 

w. 

l 

1 

1 

-0.57735 

02691 

89626 

1.00000  00000 

00000 

2 

0.57735 

02691 

89626 

1.00000  00000 

00000 
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Node  3 
at  Infinity 


(a)  Global  Representation 
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(b)  Local  Representation 


Fig.  A.l.  One-Dimensional  Infinite  Finite  Element 
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